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ABSTRACT 


Piezoelectric  transducers  are  used  for  converting  between 
electrical  and  acoustic  energy  in  such  diversified  fields  as  medical 
imaging,  underwater  acoustic  propagation,  and  non-destructive  material 
testing.  This  type  of  transducer  when  used  as  a  receiver  is  sensitive 
to  the  phase  of  the  received  pressure  wave.  The  effect  of  phase  can¬ 
cellation  on  the  measurements  of  attenuation  coefficients  of  hetero¬ 
geneous  materials  with  piezoelectric  transducers  has  been  studied 
extensively.  However,  the  fact  that  the  phase  cancellation  effect 
may  also  produce  considerable  error  in  scattering  measurements  has 
only  recently  been  recognized. 

In  this  thesis,  a  mathematical  model  for  interpretation  of  error 
arising  from  the  phase  cancellation  effect  on  backscattered  waves  is 
developed.  Computer  simulation  has  been  performed  to  elucidate  the 
influence  of  phase  cancellation  on  the  average  received  pressure  in 
terms  of  transducer  aperture  size,  number  of  scatterers,  target  range, 
and  frequency  for  five  different  scattering  configurations;  a  single 
point  scatterer;  linear,  rectangular  and  random  arrays;  and  random 
volume  distribution  of  point  scatterers. 

Results  of  the  computer  simulation  demonstrate  that  judicious 
choice  of  transducer  aperture,  target  range,  and  operating  frequency 
is  necessary  to  minimize  artifact  induced  by  the  phase  cancellation 
effect.  Based  on  a  maximum  10%  error,  aperture  values  of  .75  cm  or 
less,  farfield  target  range  and  frequency  range  of  1.00  to  15.00  mega¬ 
hertz  should  be  employed  when  implementing  phase-sensitive  transducers 


as  receivers. 
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CHAPTER  I 


INTRODUCTION 


Sound  is  a  mechanical  wave  which  can  be  produced  by  an  object 
vibrating  in  a  medium,  such  as  air,  water,  metal,  or  tissue.  Sound 
waves  have  been  found  extremely  useful  in  such  diversified  fields  as 
medical  imaging,  non-destructive  material  testing  and  underwater  war¬ 
fare.  Ultrasound  refers  to  sound  waves  above  the  frequency  of  human 
hearing  (approximately  20  kilohertz).  It  is  sensitive  to  density  and 
structural  variations  in  biological  tissue  and  materials.  In  the  case 
of  ultrasonic  medical  imaging,  in  vivo  measurements  are  possible  with¬ 
out  harmful  side  effects  of  ionizing  radiation,  such  as  X-rays.  Simi¬ 
larly,  ultrasonic  material  evaluation  for  flaw  detection  and  fatigue 
can  be  performed  without  actual  destruction  of  the  material  during 
the  testing  procedure.  Fundamentally,  scattering  and  reflection  of 
sound  waves  occurs  if  there  is  a  mismatch  of  acoustic  propagation  prop¬ 
erties  such  as  velocity,  density,  and  compressibility  between  two  media. 
The  reflected  or  scattered  sound  waves  are  detected  and  electronically 
processed  to  provide  information  regarding  the  dimensions  and  locations 
of  nonuniformities  within  the  interrogated  material.  Overall,  the 
success  of  utilizing  sound  waves  as  a  medium  for  investigation  purposes 
lies  within  the  non-destructive  nature  of  sound  and  its  sensitivity 
to  changes  in  elastic  properties  of  the  material  on  a  level  significant 
to  imaging  requirements.  In  some  applications,  sound  is  generated 
by  electrically  exciting  a  piezoelectric  crystal  such  as  quartz  or 
Lead  Zirconate  Titanate  (PZT). 


The  piezoelectric  transducer,  when  used  as  a  receiver  of  sound 
waves,  however,  is  sensitive  to  the  phase  of  the  impinging  sound  wave. 
This  phase  sensitive  nature  of  the  piezoelectric  transducer  may  give 
rise  to  serious  consequences  in  situations  involving  the  determinations 
of  acoustic  propagation  parameters  such  as  attenuation  and  scattering 
of  inhomogeneous  materials.  This  is  basically  due  to  the  non-uniform¬ 
ness  of  the  wavefront  of  the  scattered  wave  or  waves  which  have  tra¬ 
versed  through  heterogeneous  materials  arriving  at  the  face  of  a  finite- 
aperture  transducer.  This  fact  implies  that  the  phase  differences 
among  wavelets  arriving  at  different  areas  of  the  transducers  surface 
may  be  significant  enough  that  signal  cancellation  results. 

The  artifact  caused  by  the  phase-sensitive  nature  of  the  piezo¬ 
electric  receivers  which  is  termed  the  phase  cancellation  effect 
in  the  attenuation  measurements  of  an  inhomogeneous  medium  was  first 
addressed  by  Marcus  and  Carstensen  (1975).  Their  work  indicated  a 
good  correlation  between  the  relative  degree  of  sample  inhomogeneitv 
and  the  magnitude  of  the  absorption  error,  as  measured  by  a  conven¬ 
tional  phase-sensitive  receiver.  Marcus  and  Carstensen  experimentally 
compared  a  radiation  force  known  as  an  acoustic-electric  receiver  (AET) 
with  a  conventional  piezoelectric  receiver  (PET)  for  the  measurement  of 
the  absorption  coefficient  of  various  homogeneous  and  inhomogeneous 
materials.  For  2°j  agar  samples,  an  acoustically  uniform  material, 
the  attenuation  measured  by  both  types  of  transducers  proved  similar. 
However,  for  beef  muscle,  an  acoustically  inhomogeneous  tissue,  attenu¬ 
ation  values  differed.  For  measurements  at  2.00  megahertz,  the  piezo¬ 
electric  receiver  had  absorption  coefficients  ranging  from  0.21  to 
1.15  nepers/'cm,  whereas  a  radiation  force  receiver  measured  0.21  to 


0.24  nepers/cm.  The  wide  range  of  absorption  values  measured  by  a 
piezoelectric  transducer  was  attributed  to  the  phase  cancellation  arti¬ 
fact  on  the  surface  of  the  receiving  piezoelectric  element.  They, 
therefore,  suggested  the  use  of  a  phase-insensitive  acoustic-electric 
receiver  in  this  type  of  measurement. 

Further  work  by  Busse  and  Miller  (1976)  substantiated  the  exist¬ 
ence  of  phase  cancellation  effects,  noting  the  influence  of  transducer 
aperture  on  absorption  measurements  (indicating  that  a  reduction  of 
phase  cancellation  error  is  possible  with  a  decrease  in  transducer 
aperture)  and  experimentally  compared  the  PZT  and  AET.  The  AET,  util¬ 
izing  a  cadmium  sulfide  crystal,  virtually  eliminated  the  phase  cancel¬ 
lation  effect  arising  from  structurally  inhomogeneous  tissue  specimens. 
This  transducer  was  sensitive  to  the  total  power  of  the  incident  acoustic 
wave,  not  the  impinging  phasefront.  Heyman  (1979)  investigated  the 
effects  of  phase  cancellation  on  inhomogeneous  material  characteriza¬ 
tion  (largely  anisotropic  stressed  metal  samples)  by  PZT  and  AET  and 
found  a  notable  difference  in  the  attenuation  measurements.  Further¬ 
more,  materials  with  irregular  flatness  and  parallelism  also  influence 
the  degree  of  phase  cancellation. 

Reid,  Shung,  and  Kak  (1979)  have  expanded  the  concept  of  phase 
cancellation  to  include  scattered  or  reflected  waves.  Data  for  scattered 
waves "showed  that  the  scattered  strength  per  unit  volume  of  scatterers 
of  dilute  polysterene  spheres,  suspended  in  mixed  solution  of  water 
and  glycerine  and  measured  with  a  piezoelectric  transducer  of  1.00  cm 
in  diameter  located  16.00  cm  from  the  scatterers,  is  1.7  dB  lower  than 
that  obtained  with  a  piezoelectric  transducer  of  smaller  diameter 
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In  this  thesis  errors  which  may  be  introduced  due  to  the  phase 
cancellation  artifact  in  backscattering  measurements  are  examined  based 
on  computer  simulation  of  the  experimental  process.  Parameters  of 
interest  include  transducer  aperture,  number  of  scatterers,  distance 
between  the  transducer  and  scattering  medium,  and  frequency.  The 
pressure  received  by  the  transducer  is  compared  to  a  small  microprobe 
with  a  diameter  of  1.0  mm  which  can  be  approximated  as  a  point  receiver 
in  the  lower  megahertz  range.  This  investigation  also  includes  varied 
scattering  arrangements.  The  first  and  simplest  case  is  a  single  point 
scatterer  located  at  the  center  of  the  ultrasonic  beam  profile.  Thio 
arrangement  is  useful  for  two  reasons:  the  effect  of  path  length  dif¬ 
ferences  on  the  phase  of  the  pressure  wave  between  the  center  and  edges 
of  the  transducer  are  readily  apparent  and  computations  are  greatly 
simplified.  The  second  scattering  arrangement  is  a  linear  array,  or 
a  single  line  of  point  scatterers  extending  from  one  edge  of  the  trans¬ 
ducer  to  the  opposite  edge.  Rectangular  and  random  two-dimensional 
arrangements  of  point  scatterers  distributed  over  a  plane  parallel 
to  the  transducer  face  are  other  configurations  modelled.  These  two 
configurations,  along  with  a  random  volume  arrangement,  more  closely 
approximate  the  experimental  acoustic  backscattering  arrangements. 
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CHAPTER  II 


FORMULATION  OF  THE  PROBLEM 


To  facilitate  computation,  the  transducer  aperture  is  divided 
into  N  elements,  which  are  small  enough  so  that  the  pressure  received 
by  element  n  can  be  represented  by  the  pressure  received  at  the  center 
of  the  small  element  (point  n) ,  as  shown  in  Figure  1.  Assuming  that 
the  incident  wave  is  a  monochromatic  plane  wave,  we  have: 

P^  =  incident  pressure  at  scatterer  j 

=  P(0)eikRj  [1] 

where  k  is  the  wave  number  (radians/cm),  P(0)  is  the  pressure  transmitted 
by  the  transducer  and  R^  is  the  distance  along  the  axis  perpendicular 
to  the  transducer  face  to  scatterer  j  within  the  volume. 

For  conversion  from  the  three-dimensional  case  as  presented  in 
Figure  1  to  the  two-dimensional  case,  let  R^.  =  R,  where  R  is  the  target 
range,  or  let  L,  the  thickness  of  the  scattering  volume,  approach  zero. 
The  scatterers  would  then  be  placed  within  the  plane  rather  than  in 
the  scattering  volume.  Introducing  the  condition  R»D2/\,  where  D 
is  the  diameter  of  the  scatterers  and  A  is  the  wavelength,  the  scattered 
wave,  Pnj ,  at  point  n  on  the  transducer  face  due  to  a  scatterer  j  within 
the  scattering  volume,  can  be  approximated  by  the  following  equation: 


(o , i)  jw  ) 

P  .  =  P.  -J -  e  lrnjJ 

nj  l  r  .  J 

n) 


where  Sj(o,x)  is  the  scattering  amplitude  function  of  the  jth  scatterer 
and  o,£  are  unit  vectors  representing  observation  and  incident  directions. 
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respectively .  Substituting  from  equation  [1]  into  equation  [2] 
yields : 


P  . 
nj 


P(0)  S.  (6,1) 

r  . 
nj 


ik(R.+r  .) 
e  j  nj 1 


[3] 


Similarly  the  scattered  wave  at  the  center  of  the  transducer  face, 
point  o,  due  to  the  scatterer  j,  Pqj  is  given  by: 


P  . 
03 


PCO^CS.i)  ikfR+r  0 

- -H -  e  j  oj; 

roj 


[4] 


P  •  will  be  simultaneously  developed  with  P  .  for  use  as  a  reference 
oj  J  r  nj 

(this  will  be  elaborated  in  several  steps).  If  kD>>l,  we  can  further 
assume  that  scatterers  within  the  scattering  volume  are  identical  point 
scatterers  and  are  at  rest  within  the  scattering  arrangement.  For  iden¬ 
tical  point  scatterers,  S^(o,i)  is  a  constant  independent  of  6  and 
i,  therefore  Sj(o,i)=s. 


P  . 
nj 


P (0 ) s  ik(R+r  .) 
=  e  l]  nj J 

”3 


[5] 


and 


P  . 
oj 


P(0)s  eik(Rj+roj) 
roj 


[6] 


Here  a  single  scattering  process  is  assumed.  The  assumptions  generally 

hold  if  the  scatterer  concentration  is  low  and  kD^l.  Now  further 

assuming  R>>a  (a  is  diameter  of  the  transducer  aperture)  and  R>>L  (L 

is  the  thickness  of  the  scattering  volume! ,  r  .  =  r  . 

&  nj  oj 


s  R.  However, 


this  approximation  is  not  applicable  to  the  phase  term  because  very 
small  changes  in  path  lengths  produce  substantial  phase  variation. 
We  thus  obtain: 
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P  .  -  m*  eik<  w 

nj  R  J  J 


[7] 


and 


n  P(0)s  ikfR.+r  .) 
Poj  =  -  R  e  3 


[8] 


As  previously  stated,  the  intention  of  this  thesis  is  to  elucidate 
the  phase  cancellation  effect  on  the  transducer  aperture.  This  is 
accomplished  by  expressing  the  average  received  pressure  as  determined 
by  finite-aperture  transducers  in  terms  of  pressure  received  by  a  point 
receiver  located  at  point  o,  therefore  Pnj  is  normalized  with  respect 
to  PQj .  In  this  way,  the  phase  at  point  o  is  used  as  reference  and 
only  the  relative  phase  difference  between  point  n  and  point  o  is  con¬ 
sidered,  i.e., 


P  ./P  •  =  e*k^rnj  roj-* 

nr  OJ  J  J 


ei0nj 


[9] 


The  phase  value  at  the  center  becomes  0.0  degrees  as  a  result  of  this 
normalization  process  and  as  a  consequence,  the  calculated  phase  value 
over  a  certain  point  on  the  transducer  face  indicates  the  phase  of 
the  wave  at  that  point  relative  to  the  center  of  the  transducer. 
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If  there  are  M  scatterers  in  the  scattering  volume,  the  average 
received  pressure  Payg  received  by  the  transducer  is  then  given  by: 

N  M 

Pave  =  N  \  J  e  n;* 
n=l  3=1 

=  Pe10  [10] 

where  P,  8  are  the  magnitude  and  phase  of  Pave. 

The  amplitude  of  the  average  pressure  AMP  can  be  obtained  by  taking 

the  real  part  of  P  : 

r  ave 

AMP  =  Real (Pe10)  [11] 

Figure  2  depicts  the  first  four  scattering  arrangements  used  in 
the  computer  model.  The  particle  dimension  in  this  figure  is  exagger¬ 
ated  for  the  sake  of  illustrating  the  scattering  configurations.  Al¬ 
though  not  depicted  here,  the  random  volumetric  arrangement  is  similar 
to  the  random  planar  distribution,  except  that  the  scatterers  are  ran¬ 
domly  dispersed  within  a  cylinder  of  thickness  L. 

In  the  computer  model,  N,  the  number  of  elements  representing 
the  transducer  surface  area,  is  assigned  a  value  of  431. 


Figure  2.  Four  scattering  arrangements  which  are 
studied  are  graphically  illustrated. 
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CHAPTER  III 


RESULTS  AND  DISCUSSION 


The  results  and  discussion  are  presented  in  the  following  according 
to  the  five  scattering  configurations  studied. 

Single  Point  Scatterer 

The  significance  of  a  physical  phenomenon  can  usually  be  revealed 
more  readily  with  simple  examples.  Therefore,  this  study  included 
the  simplest  scattering  arrangement  possible,  depicted  in  Figure  3, 
which  shows  a  spherical  wave  originated  from  a  scatterer  located  at 
the  center  of  the  ultrasonic  beam  impinging  upon  circular  transducers 
with  aperture  sizes  of  AX  and  A2.  It  becomes  readily  apparent  from 
this  figure  that  the  phase  differences  9X  and  92  between  waves  arriving 
at  the  center  and  the  edge  of  the  transducers  depend  upon  the  aperture 
size  of  the  transducer  and  the  target  range  was  well  as  the  frequency 
of  the  wave.  These  are  the  three  parameters  which  have  been  investi¬ 
gated  for  this  case. 

Figures  4  through  9  show  the  three-dimensional  normalized  amplitude 
distribution  across  the  transducer  face  for  apertures  of  .10,  .25, 

.50,  .75,  1.00,  and  1.25  cm  at  a  frequency  of  5.00  megahertz  and  range 
of  10.00  cm.  Figure  4  establishes  that  the  pressure  amplitude  over 
a  transducer  with  an  aperture  of  1.00  mm  in  diameter  is  virtually  con¬ 
stant.  Therefore,  the  phase  cancellation  effect  for  transducers  with 
small  apertures  is  negligible.  It  will  be  seen  later  on  that  this 
statement  holds  for  all  scattering  configurations  studied.  This  is 


Figure  3.  Phase-cancellation  at  the  transducer  surfaces  of 
different  size  due  to  a  spherical  wave  generated 
by  a  single  scatterer  located  at  the  center  of 
the  incident  beam. 


Figure  5.  Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  transducer 
of  .25  cm  in  diameter  situated  at  10.00  cm  away  from  the  transducer. 


figure  6.  Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  Ml Iz  transducer 
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Figure  9.  Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  transducer 
of  1.25  cm  in  diameter  situated  at  10.00  cm  away  from  the  transducer. 
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part  of  the  reason  that  1.00  mm  aperture  probe  has  been  adopted  as 
the  reference  to  which  results  obtained  for  transducers  of  larger  size 
are  compared.  Results  for  an  aperture  of  .25  cm,  shown  in  Figure  5, 
also  indicate  an  excellent  response  with  little  error  due  to  phase 
cancellation.  Upon  increasing  the  aperture  to  .50  cm  a  slight  dis¬ 
tortion  appears  at  the  edges  of  the  transducer  resulting  in  a  moder¬ 
ate  6.88%  drop  in  the  average  received  pressure  (Figure  6).  The  maxi¬ 
mum  phase  difference  of  37.49°  occurs  near  the  edge  of  the  transducer. 
The  average  received  pressure  drops  an  overall  31.94%  as  the  transducer 
aperture  is  increased  to  .75  cm.  The  two-dimensional  display  of  the 
phase  distribution  across  this  transducer  face  is  shown  in  Figure  10. 

The  phase  value  at  the  center  of  the  transducer  is  zero  degrees  as 
discussed  in  the  section  on  the  mathematical  formulation.  To  simplify 
the  illustration  without  loss  of  meaningful  data,  only  the  integer 
values  at  each  location  were  printed.  Note  the  symmetry  of  the  phase 
distribution  across  the  transducer  face.  It  follows  that  the  ampli¬ 
tude  distribution  would  also  be  symmetrical  as  seen  in  Figure  11. 

Further  enlargement  of  the  transducer  apertures  increases  the  influence 
of  phase  cancellation  on  the  amplitude  and  phase  distributions,  average 
received  pressure,  and  phase  differences  (Figures  8  and  9).  Results 
for  an  aperture  of  1.25  cm  show  that  the  amplitude  pressure  vary  across 
the  full  range  of  1.00  to  -1.00;  similarly,  the  phase  varies  from  0.00° 
to  360.00°.  It  is  reasonable  to  predict,  just  based  on  these  results, 
that  experimental  data  measured  from  a  transducer  of  this  aperture 
(target  range  of  10.00  cm  and  frequency  of  5.00  megahertz)  would  not 
best  portray  the  acoustic  properties  of  the  material  under  investigation. 

Target  range  has  a  strong  effect  on  transducer  performance  as 
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Figure  10.  Two-dimensional  normalized  phase  distribution  on  the  face  of  a  5.00  MHz  transducer  of  .75  cm 
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well.  Intuitively,  as  R,  the  distance  between  the  transducer  and  scat- 
terers  is  increased,  the  scattered  wave  front  that  reaches  the  receiver 
should  more  approach  a  plane  wave  due  to  1/R  nature  of  scattered  spheri 
cal  waves. 

Holding  aperture  constant  at  .SO  cm  (frequency:  5.00  megahertz) 
Figures  12  through  17  demonstrate,  as  range  is  increased  from  5.00 
cm  to  20.00  cm,  the  normalized  received  amplitude  variation  on  the 
transducer  face  diminishes.  Data  also  indicate,  that  the  maximum  phase 
differences  decrease  markedly  from  74.95°  to  18.75°  for  the  previously 
mentioned  range  values.  Distances  greater  than  10.00  cm  indicate  an 
error  in  the  amplitude  distribution  of  under  5%,  which  can  be  visually 
substantiated  from  Figure  15  through  17  where  little  edge  distortion 
is  evident. 

A  summary  of  average  received  pressure  versus  target  range  at 
apertures  of  .25,  .50,  .75,  1.00,  1.25,  and  1.50  cm  is  presented  in 
Figure  18  (frequency:  5.00  megahertz).  At  an  aperture  of  1.00  cm, 
minimum  target  ranges  of  48.00  and  34.00  cm  are  necessary  for  main¬ 
taining  the  decreases  in  average  received  pressure  at  5%  and  10%  re¬ 
spectively  while  ranges  of  only  12.00  and  9.00  cm  are  required  for 
an  aperture  of  .50  cm.  Overall,  this  plot  indicates  that  for  aper¬ 
tures  above  .75  cm,  large  distances  are  required  to  reduce  the  phase 
cancellation  effect. 

The  influence  of  frequency  on  transducer  performance  is  illus¬ 
trated  in  Figure  19,  in  which  average  received  pressure  determined 
by  transducers  of  various  aperture  sizes  due  to  a  scatterer  located 
at  20.00  cm  from  the  transducer  is  plotted  versus  the  frequency  of 
the  wave.  Clearly  indicated  in  this  figure  is  that  phase  cancellation 


NORMALIZED  RECEIVED  AMPLITUDE 


Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  transducer 
of  .50  cm  in  diameter  at  a  target  range  of  10.00  cm. 
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Figure  15.  Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  transducer 
of  .50  cm  in  diameter  at  a  target  range  of  12.50  cm. 
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loinputed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  transducer 
if  .50  cm  in  diameter  at  a  target  range  of  15.00  cm. 


Average  received  pressure  measured  by  transduce 
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Figure  19.  Average  received  pressure  measured  by  transd 
.75,  1.00,  and  1.25  cm  in  diameter  due  to  a 
away  is  plotted  versus  frequency  of  the  wave 
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error  can  be  reduced  by  using  small  aperture  transducers  and  that  this 
is  especially  critical  for  frequencies  above  10.00  megahertz.  Trans¬ 
ducers  with  apertures  less  than  .50  cm  at  frequencies  under  10.00  mega¬ 
hertz  when  operated  at  this  range  suffer  a  maximum  drop  in  average 
received  pressure  of  less  than  20%. 

The  above  presented  data  for  aperture,  range,  and  frequency  demon¬ 
strate  a  clear  interrelationship  between  these  experimental  parameters 
and  the  magnitude  of  error  induced  by  the  phase  cancellation  effect. 

Linear  Scatterer  Array 

The  normalized  pressure  distribution  on  a  5.00  megahertz  transducer 

!  • 

surface  of  1.00  cm  in  diameter,  with  25 ;scatterers  arranged  in  a  linear 

format  as  indicated  in  Figure  2  located  20.00  cm  away  frqjri  the  trans- 

.  /  A 

ducer,  is  depicted  in  Figure  20.  Th'is  figure  shows  that  severe  destruc¬ 
tive  interference  due  to  phase  cancellation  occurs  near  the  edges  of 
the  transducer  perpendicular  to  the  axis  of  the  array.  Along  these 
edges  a  maximum  phase  difference  of  74.97°  was  calculated,  whereas 
at  the  opposing  edge,  where  little  error  is  evident,  a  phase  difference 
of  2.09°  was  determined.  This  amplitude  distribution  had  an  average 
received  pressure  of  .8489.  When  the  number  of  scatterers  was  increased 
to  100,  holding  all  other  parameters  constant,  only  extremely  small 
changes  in  the  three-dimensional  plot,  phase  difference  and  average 
received  pressure  occurred  (less  than  .01%).  The  three-dimensional 
plot  utilizing  100  scatterers  is  shown  in  Figure  21  and  the  associated 
average  received  pressure  is  .8485. 

Further  computation  relating  the  average  received  pressure  to 
the  mber  of  scatterers  for  5.00  megahertz  using  different  aperture 


NORMALIZED  RECEIVED  AMPLITUDE 


Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  trans 
ducer  on  the  face  of  a  5.00  MHz  transducer  of  1.00  cm  in  diameter  due  to  a  linear 
array  of  25  scatterers  located  20.00  cm  away  from  the  transducer. 


Figure  21.  Computed  normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  MHz  trails 
ducer  of  1.00  cm  in  diameter  due  to  a  linear  array  of  100  scatterers  located 
20.00  cm  away  from  the  transducer. 
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sizes  at  a  range  of  15.00  cm  (Figure  22)  reveals  that  error  arising 
from  phase  cancellation  is  independent  of  the  number  of  scatterers 
if  that  number  is  sufficiently  large.  However,  this  error  becomes 
negligibly  small  if  a  small  aperture  is  used  or  if  the  range  is  in¬ 
creased.  The  pronounced  dip  in  Figure  22  at  apertures  of  .50,  .75, 
and  1.00  cm  for  two  scatterers  indicates  the  consequences  of  severe 
phase  cancellation. 

Rectangular  and  Randomly-Distributed  Scatterers 

These  two  planar  enlargements  have  been  studied  most  extensively 
even  though  the  configuration  of  scatterers  randomly  distributed  within 
a  volume  approximates  more  closely  the  real  experimental  arrangement 
in  scattering  measurements.  This  is  because  the  computer  time  required 
for  computing  the  results  for  volume  scatterers  is  too  excessive  and 
calculations  based  on  two-dimensional  simulations  may  be  extrapolated 
to  the  three-dimensional  random  volumetric  distribution. 

Calculations  involving  the  rectangular  and  random  arrays  indicated 
that  only  a  small  number  of  scatterers  were  needed  before  a  numeric 
equilibrium  was  reached.  This  conclusion  supports  a  similar  computer 
simulation  result  presented  by  Reid,  Shung,  and  Kak  (1979).  Figure  23a 
summarizes  results  for  average  received  pressure  versus  the  number 
of  scattering  particles,  with  aperture  values  ranging  from  .25  to 
1.00  cm,  for  a  random  array  (range:  20.00  cm;  frequency:  5.00  mega¬ 
hertz).  For  all  apertures  presented,  a  minimum  of  25  scatterers  were 
needed  before  oscillations  damped  out.  This  figure  also  provides  fur¬ 
ther  evidence  that  decreasing  the  aperture  sice  increases  the  average 
received  pressure  seen  by  the  receiver,  and  indicates  apertures  equal 
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to  or  less  than  .50  cm  present  a  signal  loss  of  less  than  10 o  (average 
received  pressure  greater  than  .9000).  Figure  23b  was  generated  util- 
icing  identical  parameters  as  in  calculations  of  23a  except  that  the 
range  was  decreased  to  10.00  cm.  For  an  aperture  of  .50  cm  in  Fig¬ 
ure  23a,  the  average  received  pressure  was  approximately  .95,  whereas 
the  same  aperture  value  in  Figure  23b  yielded  an  average  received  pres¬ 
sure  of  .32.  Thus  increasing  the  distance  between  the  receiver  and 
scatterers  dramatically  improves  the  signal  strength. 

In  the  remainder  of  this  section,  the  influence  of  aperture  on 
the  phase  cancellation  effect  for  both  the  rectangular  and  random  arrays 
will  be  further  developed,  and  each  scattering  arrangement  will  be 
compared  under  identical  parameters. 

Under  the  more  realistic  conditions  of  the  rectangular  and  random 
scattering  arrangements,  increasing  aperture  distorts  the  amplitude 
and  phase  distributions  as  shown  in  Figure  24,  25,  and  26  for  the  rec¬ 
tangular  case  and  Figure  27,  28,  and  29  for  the  random  case.  For  both 
cases,  data  were  generated  at  a  frequency  of  5.00  megahertz  and  target 
range  of  10.00  cm  for  50  scatterers.  At  an  aperture  of  .50  cm.  Fig¬ 
ures  24  and  27  demonstrate  near  ideal  conditions  for  ultrasonic  measure¬ 
ment  with  an  overall  average  received  pressure  of  .9364  and  .9780  for 
rectangular  and  random  distributions,  respectively.  .An  increase  of 
aperture  to  .75  cm  introduces  a  greater  amount  of  error  as  shown  in 
Figures  25  and  28.  The  rectangular  distribution  had  a  16. 55°.s  drop 
in  average  received  pressure  while  the  random  array  distribution  had 
a  similar  drop  of  17. 61^.  The  result  of  employing  an  inappropriate  set 
of  experimental  conditions  are  Illustrated  in  Figure  26  and  29  where 
severe  amplitude  fluctuations  are  seen  especially  in  the  random  case. 
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Figure  23a.  Average  received  pressure  measured  by  transducers  of  .25,  .50, 
.75,  and  1.00  cm  in  diameter  at  5.00  megahertz  and  at  a  range 
of  20.00  cm  versus  the  number  of  scattcrers. 
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Figure  24.  Normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  megahertz  transducer 
of  .50  cm  in  diameter  due  to  50  scatterers  located  10.00  cm  away  from  the  transducer. 


Figure  27.  Normalized  pressure  amplitude  distribution  on  the  face  of  a  5.00  megahertz  transducer 
of  .50  cm  in  diameter  due  to  50  scatterers  located  10.00  cm  away  from  the  transducer. 
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The  average  received  pressure  of  each  distribution  decreased  substan¬ 
tially  (rectangular:  .7708;  random  planar:  .5255).  Up  to  an  aperture 
value  of  1.25  cm,  the  results  obtained  for  the  rectangular  and  random 
planar  distribution  were  nearly  identical. 

Increasing  the  distance  between  the  transducer  and  target 
improves  average  received  pressure.  Figure  30  depicts  average  received 
pressure  calculated  at  apertures  of  .25,  .50,  .75,  1.00,  and  1.25  cm 
with  200  scatterers  at  a  frequency  of  5.00  megahertz  utilizing  the 
rectangular  scattering  arrangement.  Apertures  greater  than  .75  cm 
need  values  of  range  greater  than  30.00  cm  before  average  received 
pressure  approaches  .9000.  .An  aperture  of  .25  cm  needs  only  a  range 
of  3.00  cm  for  a  1 0°s  signal  drop;  however,  a  .50  cm  transducer  requires 
12.50  cm.  To  obtain  an  average  received  pressure  of  approximately 
.95,  range  values  of  5.00  and  20.00  cm  are  necessary  (apertures  of 
.25  and  .50  cm  respectively). 

Increasing  frequency  has  a  pronounced  effect  on  errors  produced 
by  phase  cancellation.  This  influence  of  frequency  on  average  received 
pressure  is  shown  by  Figure  31.  The  average  received  pressure  was 
calculated  from  a  random  distribution  of  200  particles  at  a  target 
range  of  20.00  cm,  by  varying  frequency  at  apertures  of  .10,  .25,  .50, 
.75,  and  1.00  cm.  For  a  frequency  range  of  1.00  to  10.00  megahertz, 
apertures  of  .10,  .25,  and  .50  cm  provide  reasonable  performance  in 
terms  of  minimal  signal  loss  due  to  phase  cancellation.  The  data  cal¬ 
culated  for  .10  cm  transducer  aperture  further  verifies  the  near  ideal 
performance  expected  for  a  microprobe  used  as  a  receiver  and  its  role 
in  comparing  results  of  various  transducer  apertures.  When  using  aper¬ 
tures  greater  than  .25  cm  at  frequencies  above  10.00  megahertz,  the 
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Figure  30.  Average  received  pressure  measured  by  transducers  of  .25,  .50, 
.75,  1.00,  and  1.25  cm  due  to  200  scatterers  arranged  in  a 
rectangular  distribution  at  5.00  megahertz  versus  range. 


Figure  31.  Average  received  pressure  measured  by  transducers  of  .10,  .25, 
.50,  .75,  and  1.00  cm  in  diameter  due  to  200  scatterers  at  a 
range  of  20.00  cm  versus  frequency. 
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average  received  pressure  declined  rapidly  in  magnitude.  See 
Appendix  A  for  discussion  of  the  oscillations. 

Random  Volumetric  Distribution 

The  random  volumetric  scattering  arrangement  (illustrated  in 
Figure  1)  was  simulated  by  a  cylinder,  with  the  base  aligned  parallel 
to  the  transducer  face.  During  the  simulation  the  diameter  of  the 
cylinder  was  assigned  the  value  of  the  interrogating  transducer  aper¬ 
ture  for  simplification  of  geometry.  Simulation  results  depicted 
by  Curve  a  in  Figure  32  for  L  equal  to  1.00  mm,  range  equal  to 
15.00  cm,  frequency  equal  to  5.00  megahertz  and  an  aperture  of 
.25  cm,  demonstrated  that  2500  scatterers,  which  correspond  to  a 
scatterer  volume  concentration  of  approximately  500  per  cubic  milli¬ 
meter,  were  necessary  before  an  average  received  pressure  of  .9  was 
achieved.  The  identical  computation  was  performed  for  a  cylinder 
depth  of  2.00  mm  (Curve  b  in  Figure  32).  A  depth  value  of  1.00  mm 
resulted  in  a  more  rapid  convergence  to  a  low  signal  loss  condition 
than  a  depth  of  2.00  mm.  It  became  apparent  that  a  greater  number 
of  particles  is  required  to  obtain  the  same  accuracy  even  for  cases 
of  extended  range  and  small  aperture.  However,  this  requirement 
should  not  be  of  much  concern  in  biomedical  ultrasound  because  it 
is  generally  satisfied  by  biological  specimens.  For  example, 
erythrocyte  concentration  in  normal  blood  approaches  five  million 
per  cubic  millimeters. 
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CHAPTER  IV 

SUMMARY  .AND  CONCLUSION'S 

Piezoelectric  transducers  are  sensitive  to  the  phase  of  the 
incident  pressure.  Substantial  error  may  result  if  such  a  phase-sen¬ 
sitive  transducer  is  used  when  measuring  the  acoustic  properties,  such 
as  attenuation  and  scattering,  of  inhomogeneous  materials  which  dis¬ 
tort  the  shape  of  the  wave  front.  This  is  due  to  the  phase  cancella¬ 
tion  effect.  A  mathematical  model  for  the  interpretation  of  error 
arising  from  the  phase  cancellation  effect  on  the  measurement  of  back- 
scattered  waves  has  been  developed.  A  computer  program  was  implemented 
on  a  DEC- 10  to  simulate  the  ultrasonic  wave  interaction  between  a 
piezoelectric  transducer  and  an  acoustic  scattering  environment.  The 
influence  of  the  phase  cancellation  effect  on  amplitude  and  phase  dis¬ 
tributions  across  the  surface  area  of  the  ultrasonic  receiver  were 
quantified  in  terms  of  aperture  size,  frequency,  target  range,  and 
number  of  scatterers.  This  investigation  included  varied  scattering 
arrangements:  a  single  point  scatterer;  linear,  rectangular,  and  ran¬ 
dom  arrays;  and  random  volumetric  distribution. 

To  minimize  the  influence  of  phase  cancellation  on  ultrasonic 
measurements,  proper  experimental  parameters  such  as  aperture,  fre¬ 
quency,  and  range  may  be  judiciously  selected.  The  results  from  this 
investigation  arc  summarized  in  Figures  55,  54,  and  55  where  average 
received  pressure  is  plotted  versus  ratios  of  R/D,  R/F,  and  R\/D  (\ 
is  wavelength,  R  is  the  range,  and  F  is  frequency)  in  terms  of  data 
obtained  for  the  random  array  case.  Although  R/D  and  R/F  are  ulti- 
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mately  included  within  the  plot  of  average  received  pressure  versus 
RX/D,  they  are  useful  if  the  researcher  is  constrained  to  a  particular 
frequency  or  aperture  in  the  experimental  apparatus. 

Results  shown  in  Figure  33  indicate  the  necessity  of  high  R/D 
values  for  minimization  of  error  due  to  phase  cancellation.  These 
results  are  tabulated  from  data  generated  at  5.00  megahertz.  For 
other  frequencies,  separate  R/D  values  must  be  calculated.  R/D  values 
of  40.00  and  25.00  are  required  for  10%  and  5%  drops  in  average  re¬ 
ceived  pressure,  respectively.  The  units  of  range  and  transducer 
diameter  are  centimeters. 

Figure  34  represents  a  graphical  summary  of  average  received 
pressure  versus  R/F  for  apertures  of  .25,  .50,  .75,  1.00,  1.25,  and 
1.50  cm  with  the  line  demarking  an  aperture  of  .25  labelled  a  and 
continuing  through  f  for  each  respective  aperture  value.  R,  the  tar¬ 
get  range,  is  measured  in  centimeters  and  F  is  in  megahertz;  there¬ 
fore,  the  units  of  R/F  are  cm-sec/cycle.  Assuming  a  10%  drop  in 
average  received  pressure  (.9000)  and  choosing  a  diameter,  the  reader 
can  extrapolate  to  the  x-axis  to  determine  the  minimum  R/F  value. 

For  example,  selecting  a  diameter  of  1.00  cm,  a  R/F  of  no  less  than 
7.00  is  required.  For  a  larger  diameter,  a  greater  R/F  ratio  is 
necessary  to  maintain  an  adequate  average  received  pressure.  Further¬ 
more,  if  the  experiment  can  only  tolerate  a  five  percent  drop  in 
average  received  pressure  (.9500),  then  much  higher  values  of  R/F 
are  needed.  It  is  important  to  note  that  for  the  range  of  frequencies 
most  often  encountered  in  medical  imaging  (1.00  to  15.00  megahertz), 
R/F  values  greater  than  15.00  are  desirable  assuming  a  10%  pressure 
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Figure  35.  Average  received  pressure  versus  RA/D  for  the  random  planar 
distribution. 
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The  influence  of  aperture,  range,  and  frequency  on  the  phase 
cancellation  effect  is  graphically  summarized  in  Figure  35  in  which 
range  was  measured  in  centimeters,  transducer  diameter  in  centimeters 
and  wavelength  in  millimeters.  Minimum  RA/D  values  of  25.00  and  15.00 
are  necessary  to  achieve  an  experimental  accuracy  of  5%  and  10%  re¬ 
spectively.  The  units  of  RA/D  is  in  millimeters. 

In  conclusion,  the  body  of  data  provided  in  this  thesis  should 
be  very  useful  to  establishing  criteria  which  enable  investigators 
engaged  in  research  involving  the  measurements  of  acoustic  parameters, 
such  as  attenuation  and  scattering  of  heterogeneous  materials,  to 
choose  suitable  experimental  conditions  for  their  measurements  when 
phase-sensitive  devices  are  employed. 
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APPENDIX  A 


HISTOGRAM  OF  AVERAGE  RECEIVED  PRESSURE 


Investigation  into  the  variability  over  a  series  of  trials  of 
average  received  pressure  calculated  from  a  single  point  scatterer  ran¬ 
domly  placed  on  a  plane  parallel  to  the  transducer  face  provides  fur¬ 
ther  evidence  of  the  influence  of  increased  aperture  on  the  phase 
cancellation  effect.  The  data  was  calculated  by  the  FORTRAN  program 
ARPVAR.  Results  were  plotted  as  a  histogram  for  each  modelling  condi¬ 
tion. 

For  25  trials,  the  average  received  pressure  was  calculated  on 
transducer  apertures  of  .50,  .75,  1.00,  and  1.25  cm  (Figures  36,  37, 

38,  and  39)  at  a  range  of  7.50  cm  and  a  frequency  of  5.00  megahertz. 

For  an  aperture  of  .50  cm  received  pressures  for  all  25  trials  fell 
within  the  .90  to  1.00  interval.  The  mean  average  received  pressure 
was  .9798  and  had  a  variance  of  .0006.  Increasing  the  aperture  to 
.75  cm  resulted  in  a  decrease  of  the  mean  average  received  pressure 
to  .9008  and  increased  sample  variance  to  .00138.  Referring  to  Fig¬ 
ure  38,  the  variability  of  average  received  pressure  outcomes  in¬ 
creased  with  increasing  aperture.  For  1.00  cm  aperture,  the  sample 
distribution  had  a  mean  of  .7148  and  variance  of  .00961.  Concluding 
with  Figure  39,  an  aperture  of  1.25  cm  shows  the  large  fluctuations 
in  average  received  pressure  due  to  phase  cancellation  calculated  for 
a  single  point  scatterer.  This  outcome  is  expected  and  emphasises 
that  correct  aperture  size  is  extremely  important  when  utilising  phase- 
sensitive  piezoelectric  transducers. 


The  variability  of  average  received  pressure  was  investigated 
at  a  range  of  20.00  cm,  a  frequency  of  25.00  megahertz  and  200  scat- 
terers  for  apertures  of  .50,  .75,  and  1.00  cm  (Figures  40,  41,  and 
42) .  Twenty-five  trials  were  performed  at  each  of  these  aperture 
values.  These  parameters  were  chosen  in  order  to  compare  data  sum¬ 
marized  in  Figure  51.  In  Figure  31,  each  point  was  obtained  from 
only  one  trial.  As  evident  from  Figures  41  and  42,  there  is  consider¬ 
able  variation  from  trial  to  trial  for  apertures  of  .75  and  1.00  cm 
and  results  of  one  trial  are  not  representative  of  the  typical  out¬ 
come  under  these  conditions.  For  example,  the  results  in  Figure  42 
demonstrate  a  large  variance  in  average  received  pressure  when  using 
a  1.00  cm  aperture  transducer.  The  sample  mean  of  this  distribution 
is  .2525  and  has  a  variance  of  .0463.  The  consequence  is  that  re¬ 
sults  obtained  from  one  trial  may  differ  substantially  from  the 
other.  Therefore,  data  for  higher  frequencies,  e.g.,  data  above 
13.00  megahertz  at  an  aperture  of  1.00  cm  and  above  22.00  mega¬ 
hertz  at  an  aperture  of  .75  cm  would  be  expected  to  show  considerable 
variation  for  different  arrangements  of  200  scatterers.  This  result 
explains  the  large  fluctuations  seen  in  Figure  31  at  higher  frequen¬ 
cies  for  large  apertures. 


Figure  36.  Average  received  pressure  versus  the  number  of  occurrences 
for  a  single  point  scatterer  at  a  range  of  7.50  cm  and 
aperture  of  .50  cm  over  25  trials. 


Figure  38.  Average  received  pressure  versus  the  number  of  occurrences 
for  a  single  point  scatterer  at  a  range  of  7.50  cm  and 


Average  received  pressure  versus  the  number  of  occurrences 
for  a  single  point  scatterer  at  a  range  7.50  cm  and 
aperture  of  1.25  cm  over  25  trials. 
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APPENDIX  « 


HOW  Tr  'JSE  THE  COMPUTER  PROGRAMS 


The  computer  programs  used  in  this  simulation  were  written  to 
allow  maximum  flexibility  and  accurate  presentation  of  data  generated 
in  the  modelling  procedure.  The  author  made  extensive  use  of  modular 
program  construction  so  that  other  workers  can  clearly  understand 
program  execution  and  implement  source  code  modification.  The  oper¬ 
ating  requirements  and  logic  flow  of  each  computer  program  are 
presented. 


Phase 

A  FORTRAN  program  PHASE  was  developed  for  use  in  the  computer 
simulation.  This  program  allows  the  user  to  select  the  following 
scattering  arrangements:  a  single  point  source,  linear  array, 
rectangular  array,  random  planar  array,  and  random  volume.  The  user 
also  specifies  aperture,  number  of  scatterers,  range  and  frequency 
for  calculations  of  received  amplitude  and  phase  distributions  across 
the  surface  area  of  the  ultrasonic  transducer.  Once  the  program 
PHASE  has  been  loaded,  the  program  will  prompt  the  user  for  input 
parameters  as  shown  by  Figure  43.  Initial  calculations  are  displayed 
on  the  computer  terminal  for  verification  by  the  user.  The  result¬ 
ing  data  calculated  can  be  printed  on  paper  and/or  written  to  a  disk 
file  for  plotting  by  PLT3D,  CPLOT,  or  other  available  x-y  plotting 


routines. 
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p 1 !  a  a -  c  a  c  i :  l  la t ion  p  k  o  g  k  a  :•*. 


SCAT  T  K  R  E  R  AR  R  ANG  EM  E  NT(  L,  0,  K,  V  ,  A)  : 

SPECIFY  NUMBER  CF  SCATTERING  PARTICLES: 

ENTER  DISTANCE  FROM  TRANSDUCER  TO  PARTICLES: 


TRANSDUCER  APERTURE : 

PLEASE  ENTER  OPERATING  FR  EfJ  UE  MC  Y  (  N  EC  A  ME  RTF.  )  : 


VARY  RANGE,  SCATTERF.RS,  APERTURE,  F  REOUC  NCY  (  P.  ,  N ,  A  ,  F  ) 
ENTER  STOP  V A L U E ( REAL) : 

INCREMENT : 

DO  YOU  WANT  PLOTTINC  DATA  WRITTEN  TO  DISK(Y/N): 
PRINTED  DATA  WRITTEN  TO  DISE(Y/N): 

IS  ABOVE  INPUT  DATA  COR RECT ( Y / N ) : 


Figure  43.  Sample  input  to  program  PRASE, 


w 
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eft********************************************************************** 

c 

C  Bioengineering:  Phase-Cancellation  Studies 
C 

C  SCATTERERING  ARRANGEMENTS: 

C  1)  LINEAR 

C  2)  RECTANGULAR  ARRAY 

C  3)  RANDOM  ARRAY 

C  4)  RANDOM  VOLUME 

C  5)  RANDOM  ARRAY  CONSTANT  SCATTERING  APERTURE 

C 

Q* ********* ****** **★ *** ************************* ** * ********************* 

COMPLEX  PRS(41,4l) 

COMMON  / B/RANGE , NUM , APERT , FREQ , WAVEL, XK, PARRAD , NCOUNT, PCT, PSF , 

1  DIV , XNEAR, IF LAG , IFLAG1 , IFLAG2 , IFLAG3 , VELCTY , XKA, RFREQ( 100) , 

2  TRANR, AMP(41 ,41 ) ,THETA( 41,41) , ASF, RPHASE( 100) ,RREAL( 100) , 

3  SCAT, CENTER, M,AMPMIN, AMPMAX,THEMIN,THEMAX,PIC(41 , 41 ) ,DIM, 

4  RNUM( 100) .RPRSC100) ,RMAG(100) ,RIMAG(100) .DAMP, DTHETA, RPCT( 100) , 

5  TEMPRS , ANS , RAMP (100), STOP , XINC , SCTANS , XNUM , RVOL( 1 00 ) , ARPVOL , PRS 

C  MAIN  PROGRAM 

TYPE  1000 

1000  FORMAT(26X,' BIOENGINEERING  PROGRAM',//) 

CALL  ASK 
1005  CALL  INIT 

CALL  CRUNCH 
CALL  WARN 
CALL  ARYPIC 
CALL  DECIDE 
CALL  VERIFY 
CALL  NORMAL 
CALL  STAT 
CALL  AVGPRS 
CALL  FILINC 
CALL  INCVAR 

IF( IFLAG2 . EQ.O)  GOTO  1005 
2050  TYPE  2055 

2055  FORMAT('ODO  YOU  WISH  TO  CONTINUE(Y/N) ?  ',$) 

CALL  BELL 

ACCEPT  2060, ANSWER 
2060  FORMAT(Al) 

IF( ANSWER. NE. 'Y' )  GOTO  2070 
CALL  ASK 
GOTO  1005 

2070  CALL  PARCNT 
CALL  SAVDAT 
CALL  VOLDAT 
STOP 
END 


o  o 


C************** ****************************************** ******** 


PRTDAT :  PRINT  RESULTS  OF  AMPLITUDE  AND  PHASE  CALCULATIONS 
C 

C*  ************************************************* ************** 

SUBROUTINE  PRTDAT( FNAME3 , FNAME4 ) 

DOUBLE  PRECISION  FNAME3 , FNAME4 
INTEGER  ITHETA( 41,41) 

CALL  OPEN( 21 , 'FILE' , FNAME3) 

CALL  OPEN(22, 'FILE' ,FNAME4) 

WRITE(21 ,2076) 

2076  FORMATC 5 3X,' AMPLITUDE  DATA') 

WRITE( 22 , 2077 ) 

2077  FORMAT(58X,' PHASE  DATA') 

WRITEC 2 1,2078) 

2078  FORMAT(53X,14( '-'),//) 

WRITE( 22 , 2079 ) 

2079  FORMAT(58X,10( '-'),//) 

DO  2090  I  =1  ,M 
DO  2091  J  =1  ,M 

2091  ITHETA( I , J )  =  INT(THETA( I, J) ) 

2090  CONTINUE 

WRITE(  21 , 2085)  (  (  AMP(  I,  J)  ,  J=*9 , 33)  ,  1=9 ,33) 

WRITEC  22 , 2086 )  ( ( ITHETA(  I ,  J)  ,  J=*9 , 33 ) ,  1=9 , 33 ) 

2085  FORMAT(25(F5.2)) 

2086  FORMAT( 25(14, IX)) 

WRITEC 2 1,2081) 

WRITE(22,2081) 

2081  FORMATC'  ',///) 

IF(SCTANS.EQ.'O')  WRITE(21 , 2093)  PCT, FREQ 
IFCSCTANS.EQ.'O' )  GOTO  2094 

2093  FORMAT( IX, 'NUMBER:' ,F8.1,7X, 'FREQUENCY:' ,F5. 2,'  MHZ') 

WRITE( 2 1,2082)  NUM, FREQ 

WRITE( 22, 2082)  NUM, FREQ 

2082  FORMATClX, 'NUMBER:', 15, 10X, 'FREQUENCY:', F5 . 2,'  MHZ') 

2094  WRITEC 2 1,2083)  RANGE 
WRITEC 2 2, 2083)  RANGE 

2083  FORMATC IX, 'DISTANCE  FROM  TRANSDUCER  TO  SCATTER: ' ,F7 .2, '  CM') 
WRITEC 2 1,2084)  APERT 

WRITE (2 2, 2084)  APERT 

2084  FORMATC IX, 'TRANSDUCER  APERATURE: ' ,F6.3, '  CM') 

WRITEC 2 1,2088)  DAMP 

WRITEC 22, 2088)  DTHETA 
2088  FORMATC IX, ' DIFFERENCE: ' , FI  1 . 5) 

WRITEC 2 1,2092)  TEMPRS 

2092  FORMATC IX, 'AVERAGE  RECEIVED  PRESSURE: ' ,F10.5) 

CALL  CLOSEC21) 

CALL  CLOSEC22) 

RETURN 

END 


non 


Q*  ji  **************************  *  *****  Irk  'k'k'k'k  *  irk  rt  *  x  *  k  *  x  x  *  ***  *  tit  w*  x  ***  ** 

PRTDAT :  PRINT  RESULTS  OF  AMPLITUDE  AND  PHASE  CALCULATIONS 

Q**************************************************************** 

SUBROUTINE  PRTDAT( FNAME3 , FNAME4) 

DOUBLE  PRECISION  FNAME3 , FNAME4 
INTEGER  ITHETA( 41,41) 

CALL  OPEN(21, 'FILE' ,FNAME3) 

CALL  OPEN(22, ' FILE' , FNAME4) 

WRITE(21 , 2076) 

2076  FORMAT(53X, 'AMPLITUDE  DATA') 

WRITE( 22 , 2077 ) 

2077  FORMAT(58X,' PHASE  DATA') 

WRITE( 2 1,2078) 

2078  FORMAT(53X,14( '-'),//) 

WRITE( 22 , 2079 ) 

2079  FORMAT(58X,10( '-'),//) 

DO  2090  I  =1 ,M 

DO  2091  J  =1 ,M 

2091  ITHETA(I.J)  =  INT(THETA( I, J) ) 

2090  CONTINUE 

WRITEC21 , 2085)(( AMP( I, J) , J=9 , 33) ,1=9 ,33) 

WRITE( 2  2 , 2086 ) ( ( ITHETA( I, J) ,J*9,33) ,1=9,33) 

2085  FORMAT(25(F5.2)) 

2086  FORMAT( 25(14, IX)) 

WRITE( 2 1,2081) 

WRITE(22,2081) 

2081  FORMATC'  ',///) 

IF( SCTANS.EQ. ' O' )  WRITE(21 , 2093)  PCT.FREQ 
IF(SCTANS.EQ.'0' )  GOTO  2094 

2093  FORMATC  IX, 'NUMBER:'  ,F8.1,7X,’  'FREQUENCY:'  ,F5. 2,'  MHZ') 
WRITE( 2 1,2082)  NUM.FREQ 

WRITE(22,2082)  NUM.FREQ 

2082  FORMATC 1 X, ' NUMBER :', 15 ,1 OX, 'FREQUENCY :',F5. 2,'  MHZ') 

2094  WRITEC 2 1,2083)  RANGE 
WRITE(22, 2083)  RANGE 

2083  FORMATC IX, 'DISTANCE  FROM  TRANSDUCER  TO  SCATTER: ' ,F7. 2) 
WRITEC 2 1,2084)  APERT 

WRITEC 22, 2084)  APERT 

2084  FORMATC IX, 'TRANSDUCER  APERATURE: ' , F6 .3 , '  CM') 

WRITEC 2 1,2088)  DAMP 

WRITEC 22, 2088)  DTHETA 
2088  FORMATClX, 'DIFFERENCE:' ,F11. 5) 

WRITEC 2 1,2092)  TEMPRS 

2092  FORMATC IX, 'AVERAGE  RECEIVED  PRESSURE: ' ,F10. 5) 

CALL  CLOSEC21) 

CALL  CLOSEC22) 

RETURN 

END 
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C******************************************************************** 

c 

C  AVGPRS:  CALCULATE  AVERAGE  RECEIVED  PRESSURE  FOR  THE  TRANSDUCER 
C 

c* ************************************ ************************* ****** 
SUBROUTINE  AVGPRS 
K  =  0 

TEMPRS  =0.0 

DO  3005  I  =  1,M 

DO  3000  J  =  1 ,M 

IF(PIC( I, J) .EQ. '-' )  GOTO  3000 

TEMPRS  =  TEMPRS  +  AMP(I,J) 

K.  =  K  +  1 
3000  CONTINUE 
3005  CONTINUE 

TEMPRS  =  TEMPRS/ FLOAT(K) 

RPRS(NCOUNT)  =  TEMPRS 

AREA  =  3.1415927  *  (TRANR**2) 

APRESS  =  (TEMPRS/ AREA) 

TYPE  30 10, APRESS .TEMPRS 

3010  FORMAT( 'OAVGPRS/AREA  =  ' ,F11 .5, '  AVGPRS  =' ,F11 .5) 
IF(SCTANS.NE.'V')  GOTO  3025 
VOLUME  =  3.1415926  *  (TRANR**2)  *  .1 
ARPVOL  =  TEMPRS/VOLUME 
RVOL(NCOUNT)  =  ARPVOL 
3025  RETURN 
END 

C******************************************************************** 

C 

C  CRUNCH:  CALCULATE  SCALE  FACTORS 
C 

C*********************** ************************* ******************** 

SUBROUTINE  CRUNCH 
NNUM  =  NUM 

IF(NUM.EQ.l)  NNUM  =  NNUM  +  1 
XNUM  =  FLOAT(NNUM) 

PSF  =  ( APERT/ ( XNUM- 1 .0)) 

CENTER  =  INT(M/2 .0)  +  1.0 
TRANR  =  (APERT/ 2.0) 

ASF  =  (APERT/ ( DIM- 1 .00)) 

PARRAD  =  6.0 

RETURN 

END 


o  n 
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c**************************************************************** 

C 

INFO:  WRITE  INFORMATION  TO  DISK  IN  FILE  INFXXX.DAT 


c**************************************************************** 


SUBROUTINE  INFO( FNAME1 , FNAME2) 

DOUBLE  PRECISION  FNAME1 , FNAME2 
CALL  OPEN( 21 , ' FILE' , FNAME2) 

WRITE! 2 1,4000) 

4000  FORMAT!  26X,' INFORMATION  FILE') 

WRITE( 2 1,4005) 

4005  FORMAT(26X,16( '-'),///) 

WRITE( 2 1,4010)  FNAME1 

4010  FORMAT(5X,' PLOTTING  DATA  FILE  NAME:  ' ,A10,/) 

WRITE(21 , 4015)  NUM 

4015  FORMATC5X,' NUMBER  OF  SCATTERING  PARTICLES 15 , /) 

WRITE( 2 1,4020)  RANGE 

4020  FORMATC5X,' DISTANCE  TRANSDUCER-  SCATTERERS: ' , F7 .2, IX, ' CM' ,/ ) 
WRITE( 2 1,4025)  APERT 

4025  FORMAT(5X, ' APERATURE  OF  TRANSDUCER: ' ,F8.3, IX, ' CM' ,/) 

WRITE( 2 1,4035)  FREQ 

4035  FORMAT! 5X, ' FREQUENCY : ' , F7 . 2 , IX, ' MEGAHERTZ ' , / ) 

WRITE( 2 1,4040)  WAVEL 

4040  FORMATC5X, ' WAVELENGTH: ' ,F7. 2, IX, 'MILLIMETERS' ,/) 

WRITE( 2 1,4045)  XK 

4045  FORMAT! 5X,' WAVENUMBER  K:' ,F10. 2,'  RAD/CM',/) 

WRITE( 2 1,4050)  PARRAD 

4050  FORMAT(5X, ' SCATTERER  RADIUS: ' ,F7 .2, IX, 'MICROMETERS' ,/) 

WRITE! 2 1,4060)  DIV 

4060  FORMAT! 5X,' ANGLE  OF  DIVERGENCE: ' ,F7 .2, IX, ' DEGREES' ,/) 

WRITE! 2 1,4065)  XNEAR 

4065  FORMAT(5X, ' NEAR-FARFIELD  BOUNDARY :', F7 . 2 , IX, ' CENTIMETERS' ,/ ) 
WRITE( 2 1,4070)  VELCTY 

4070  FORMAT(5X, 'WAVE SPEED : ' ,F8 .2, IX, 'METERS/ SEC' ,/) 

WRITE(21 ,4075)  XKA 

4075  FORMAT! 5X, ' KA: ' , F8 . 5 , /) 

WRITE( 2 1,4080)  SCAT 

4080  FORMAT(5X,' SCATTERING  PARTICLE  COEFFICIENT: ' ,F5.2,/) 

WRITE! 2 1,4090)  SCTANS 

4090  FORMAT(5X,'SCATTERERING  ARRANGEMENT:  ',A1,/) 

WRITE(21 ,4100)  TEMPRS 

4100  FORMAT(5X, 'AVERAGE  RECEIVED  PRESSURE  =',F11.5,/) 

WRITE(21 ,4105)  DAMP 

4105  FORMAT! 5X, 'AMPLITUDE  DIFFERENCE  -  TRANSDUCER  =  ' , FI  1 . 6 , / ) 

WRITE! 2 1,4 110)  DTHETA 

4110  FORMAT! 5X, 'PHASE  DIFFERENCE  ACROSS  TRANSUCER  =  ' ,F11.6,/) 

CALL  CLOSE! 21) 

RETURN 

END 


(]* ****************************************************** ********* 

c 

C  VERIFY:  DISPLAY  INFORMATION  ON  TTY  FOR  VERIFERCATION 
C 

C************************** ****** ******************************** 

SUBROUTINE  VERIFY 
TYPE  5000 

5000  FORMAT( ' 0  ') 

TYPE  5005 

5005  FORMAT(26X,' INPUT  DATA  SUMMARY' ,/) 

TYPE  5010 , NCOUNT, PCT 

5010  FORMAT('ONCOUNT:' ,16, 25X, ' NUMBER  OF  PARTICLES' , F7 . 1 ) 

TYPE  5015, RANGE, APERT 

5015  FORMAT( 'OLENGTH  FACE  TO  PARTICLES' ,F7 .2, 7X, 'APERATURE: 'F8. 3) 
TYPE  5020,WAVEL,XK 

5020  FORMAT('OWAVELENGTH:',F7. 2, 20X, 'WAVENUMBER:', F10. 2,'  RAD/CM') 

TYPE  5025, PARRAD, FREQ 

5025  FORMAT( ' OPARTICLE  RADIUS: ' ,F7 .2, 15X, 'FREQ: ' ,F7 .2, IX, 'MHERTZ' ) 
TYPE  5030,DIV,XNEAR 

5030  FORMAT( 'ODIVERGENCE: ' ,F7 .2, 20X, ' NEARFIELD  BOUNDARY :', F7 . 2 ) 
RETURN 
END 

C********************* ******************** ************************ 

c 

C  INIT:  INTIALIZE  ULTRASONIC  AND  ARRAY  PARAMETERS 
C 

C* ***************************** *********************************** 

SUBROUTINE  INIT 

M  =  41 

DIM  =  23.0 

SCAT  =  1.0 

VELCTY  =  1500.00 

WAVEL  =  (VELCTY/ ( FREQ*1 .0E03) ) 

XK  =  (2.0  *  3.14159)/ (WAVEL/ 10.0) 

XNEAR  =  (APERT**2)/ ( .40  *  WAVEL) 

DIV  =  ( ASIN( ( . 122  *  WAVEL) /APERT))  *  57.29578 

NCOUNT  =  NCOUNT  +  1 

AMPMIN  =  1000.00 

AMPMAX  =-1000.00 

THEM IN  =  1000.00 

THEMAX  =-1000.00 

RETURN 

END 


o  o  o  o 
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C* *************************************************************** 


SAVDAT:  CREATE  DATA  FILE  FOR  PLOTTING  OF  AVERAGE  AMPLITUDE 
VS.  NUMBER  OF  SCATTERERS  AND  WRITE  DATA  SUMMARY  FILE  TO  DISK 

C* ********** ******************************************** ********* 

SUBROUTINE  SAVDAT 

DOUBLE  PRECISION  FNAME5 , FNAME6 

CALL  NAME( FNAME5 , FNAME6 ) 

CALL  OPEN(21 , ' FILE' , FNAME5) 

CALL  0PEN(22,'FILE'  ,FNAME6) 

WRITE(22,8010) 

8010  FORMAT( 35X, ' DATA  SUMMARY') 

WRITE(22,8015) 

8015  F0RMAT(35X,12( '-'),////) 

WRITE(22,8016)  SCTANS 

8016  F0RMAT(6X, ' SCATTERERING  ARRANGEMENT:  ',Al,//> 

IF(ANS.NE.'F')  GOTO  8022 

WRITE(22,8021)  RANGE, APERT, NUM 

8021  FORMAT(6X,  'RANGE : ' ,F7 .2, '  CM  APERT: ' ,F6.3, '  CM  NUM:  ',15,//) 
WRITE(22,8026) 

8026  FORMAT(5X,'FRQ' ,4X,'AVG  PRS' ,4X, 'DIFF  AMP' , 4X, ' DIFF  PHASE', 4X, 

1  ' REAL' , 5X, ' IMAGINARY' , 4X, ' MAGNITUDE' ) 

GOTO  8002 

8022  WRITE(22,8020)  RANGE, APERT, FREQ 

8020  FORMAT(6X, ' RANGE: ' , F7 . 2, '  CM  APERT: ' ,F6.3, '  CM  FREQ:  ', 

1  F5.2, '  MEGAHERTZ',//) 

WRITE(22,8025) 

8025  FORMAT( 5X, ' NUM' , 4X, ' AVG  PRS' ,4X, 'DIFF  AMP' , 4X, ' DIFF  PHASE', 4X, 

1  ' REAL' , 5X, ' IMAGINARY' , 4X, ' MAGNITUDE' ) 

8002  WRITE(22,8030) 

8030  FORMAT(5X,' - ' ,4X, 7( '-' ) ,4X, 8( ) ,4X, 10( '-' ) ,4X, ' - ',5X, 

1  9('-'),4X,9 ('-'),//) 

DO  8000  I  =1 .NCOUNT 
IF(ANS.NE.'F')  GOTO  8001 
WRITE( 2 1,8005)  RFREQ( I) , RPRS( I) 

WRITE(22,8040)  RFREQ( I),RPRS(I) ,RAMP(I) ,RPHASE(I) ,RREAL(I) , 

1  RLMAG(I)  ,RMAG(I) 

GOTO  8000 

8001  WRITE(21 , 8005)  RNUM( I) , RPRS( I) 

8035  WRITE(22,8040)  RNUM( I) , RPRS( I) , RAMP( I) , RPHASE( I) , RREAL(I) , 

1  RIMAG(I) ,RMAG(I) 

8005  FORMAT(3X,F11.4,2X,F11.7) 

8040  FORMAT(3X,F6.0,2X,F9.5,3X,F7.3,7X,F7.3,4X,F7.3,5X,F7.3,5X,F7.3) 

8000  CONTINUE 

CALL  CL0SE( 21) 

CALL  CL0SE( 22) 

RETURN 

END 


n  o 


*  mm 


c****************************************************************** 

c 

C  ARYPIC:  INTIALIZE  PICTATORIAL  TRANSDUCER  REPRESENTATION 
C  ROUTINE  LOCATES  PARTICLE  SCATTERERS  WITHIN  TRANSDUCER  BEAM 
C 

c* ***************************** ************************************ 

SUBROUTINE  ARYPIC 

DO  6005  I  =  1 ,M 
DO  6000  J  =  1,M 
PIC(I,J)  = 

AMP( I, J )  =  0.0 
6000  THETA( I , J )  =  0.0 

6005  CONTINUE 
RETURN 
END 

C* ********************************** ******************************* 

C 

PLTDAT :  WRITE  PLOTTING  DATA  TO  DISK 

C* *********************************** ****************************** 


SUBROUTINE  PLTDAT(FNAMEl) 

DOUBLE  PRECISION  FNAME1 

CALL  OPEN( 21,'  FILE'  , FNAME1 ) 

WRITE( 21,7000)  M, AMPMIN, AMPMAX, THEMIN, THEMAX, SCTANS 
7000  FORMAT( 13, IX, FI 1.5, IX, FI  1.5, IX, FI  1.5, IX, FI  1.5, IX, Al) 
WRITE (2 1,7005)  NUM , TRANR, RANGE , FREQ , APERT, PCT 
7005  F0RMAT(I5,1X,F11.5,1X,F11.5,1X,F11.5,1X,F11.5,1X,F11.5) 

WRITE( 21 , 7010) ( ( AMP( I, J) , J=1 ,M) , 1=1 , M) 

WRITE(21 ,7010) ( (THETA( I, J) , J=L,M) , 1=1 ,M) 

7010  FORMAT(41 (FI 1 . 6, IX) ) 

CALL  CLOSE( 21) 

RETURN 

END 


o  o  o 
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C******************************************************************* 

ASK:  PROMPT  USER  FOR  INPUT  VARIABLES 

C******************************************************************* 


SUBROUTINE  ASK 

8000  IFLAG  =  0 
IFLAG1  =  0 

8002  TYPE  8003 

8003  FORMAT( 'OSCATTERER  ARRANGEMENT L,0,R,V, A) :  ',$) 

ACCEPT  8055, SCTANS 

TYPE  8001 

8001  FORMAT 'OSPECIFY  NUMBER  OF  SCATTERING  PARTICLES:  ',$) 

ACCEPT  8005, NUM 

8005  FORMAT  I) 

TYPE  8010 

8010  FORMAT ' OENTER  DISTANCE  FROM  TRANSDUCER  TO  PARTICLES:  ',$) 
ACCEPT  8015, RANGE 
8015  FORMAT  F) 

TYPE  8020 

8020  FORMAT 'OTRANSDUCER  APERTURE:  ',$) 

ACCEPT  80 15, APERT 
TYPE  8030 

8030  FORMAT ' OPLEASE  ENTER  OPERATING  FREQUENCY ( MEGAHERTZ ) :  ',$) 
ACCEPT  80 15, FREQ 
8045  TYPE  8050 

8050  FORMAT' ORANGE, NUMBER  SCATTERERS, APERTURE, FREQ(R, N, A, F) :  ',$) 
ACCEPT  8055, ANS 
8055  FORMAT  Al) 

IF( ANS.NE. ' R' .AND. ANS. NE. 'N' .AND. ANS. NE. ' A' .AND. ANS. NE. ' F' ) 

1  GOTO  8045 
TYPE  8070 

8070  FORMAT 'OENTER  STOP  VALUE(REAL  VALUE):  ',$) 

ACCEPT  8015, STOP 

IF( STOP. NE. FLOAT NUM))  GOTO  8074 

XINC  =  1.0 

GOTO  8079 

8074  TYPE  8075 

8075  FORMAT' 0 INCREMENT  OF  VARIABLE ( REAL  VALUE):  ',$) 

ACCEPT  80 15, XINC 

8079  TYPE  8080 

8080  FORMAT 'ODO  YOU  WANT  PLOTTING  DATA  WRITTEN  TO  DISK(Y/N)?  ',$) 
ACCEPT  8055, ANSWER 

IF( ANSWER. EQ.'N')  IFLAG  =  1 
TYPE  8085 

8085  FORMAT 'OPRINTED  DATA  WRITTEN  TO  DISK(Y/N)?  ',$) 

ACCEPT  8055, ANSWER 

IF( ANSWER. EQ.'N')  IFLAG 1  =  1 

TYPE  8090 

8090  FORMAT '0 IS  ABOVE  INPUT  DATA  CORRECT(Y/N) ?  ',$) 

ACCEPT  8055, ANSI 

IF( ANSI. EQ.'N')  GOTO  8000 

RETURN 

END 
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C*** ***************************************************** *********** 

c 

NORMAL:  THIS  SUBROUTINE  NORMALIZES  DATA  AS  OUTLINED. 

1.  DETERMINE  PHASE  FROM  COMPLEX  PRESSURE  ARRAY 
CALCULATE  AMPLITUDE  FROM  MAGNITUDE*COS( PHASE) 

LOCATE  AMPLITUDE  AND  PHASE  MAXIMA  AND  MINIMA 
WRITE  RESULTS  TO  ARRAYS 


C 
C 
C 
C 
C 
C 

***************************************************** *********** 


2. 

3. 

4. 


SUBROUTINE  NORMAL 
COMPLEX  PRS(41 , 41 ) ,Z 


Z  =  PRS(IFIX( CENTER), IFIX( CENTER)) 


DO  7005  I  =  10,32 

DO  7000  J  =  10,32 

IF(  P  IC(  I ,  J  )  .  EQ . '  - '  )  GOTO  7000 

PRS( I, J)  =  PRS(I, J)/Z 

THETA( I, J)=57 . 295780*ATAN2( AIMAG(PRS(I, J) ) , REAL(PRS( I, J) ) ) 
AMP( I, J)  =  COSD(THETA( I, J) ) 

IF(THETA(I, J) .LT.0.0)  THETA(I.J)  =  THETA( I , J )  +  360.00 
7005  CONTINUE 
7000  CONTINUE 

DO  7025  I  =1 ,M 
DO  7020  J  =1 ,M 

IF(PIC(I, J) . EQ.'-' )  GOTO  7020 
IF(AMP(I, J) .GT.AMPMAX)  AMPMAX  »  AMP(I.J) 

IF( AMP( I, J) .LT. AMPMIN)  AMPMIN  =  AMP(I,J) 

IF ( THETA( I , J ) . GT . THEMAX)  THEMAX  *  THETA( I , J ) 

IF(THETA( I, J) .LT.THEMIN)  THEMIN  =  THETA(I.J) 

7020  CONTINUE 
7025  CONTINUE 

DAMP  =  AMPMAX  -  AMPMIN 
DTHETA  =  THEMAX  -  THEMIN 
RPCT(NCOUNT)  =  PCT 
RPHASE(NCOUNT)  =  DTHETA 
RAMP(NCOUNT)  *  DAMP 
RNUM(NCOUNT)  =  FLOAT(NUM) 

RMAG(NCOUNT)  =  CABS(Z) 

RIMAG(NCOUNT)  =  AIMAG(Z) 

RREAL(NCOUNT)  =  REAL(Z) 

RFREQ(NCOUNT)  =  FREQ 

RETURN 

END 


i 


o  o  o  o 
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£********************************************************************** 

RDNAPT:  THIS  SUBROUTINE  RANDOMLY  PLACES  A  POINT  SCATTERER  WITHIN  A 
SURFACE  AREA  SPECIFIED  IN  VARIABLE  SCTDIM.  THE  USER  VARY  APERTURE 

C************* *************** ****************************************** 


SUBROUTINE  RDNAPT 
COMPLEX  C,PRS(41,41) 

PCT  =  FLOAT(NUM) 

SEED I  =  .012345 
SEED2  =  .543210 
SEED3  =  .975310 
SEED4  =  .013579 
SCTDIM  =  .75 

ASFSCT  =  SCTDIM/ (DIM  -  1.0) 

SCTRAD  =  SCTDIM/ 2.0 
DO  1010  II  =  1 ,M 
DO  1005  J1  =  1  ,M 
XII  =  FLOAT(Il) 

XJ1  -  FLOAT(Jl) 

TEMPI  =  ASF*SQRT((XI1-CENTER)**2  +  ( XJ 1-CENTER)** 2 ) 

IF( TEMPI. GT.TRANR)  GOTO  1005 
C  -  (0.0, 0.0) 

NUMBER  =  0 
PIC(I1,J1)  =  'X' 

DO  1000  12  -  1,100000 
SIGN1  =1.0 
SIGN2  =1.0 
TEMP5  =  RAN(SEED3) 

TEMP6  =  RAN(SEED4) 

IF(TEMP5.LT. .5000)  SIGN1  =  -1.0 

IF(TEMP6.LT. .5000)  SIGN2  =  -1.0 

XRAND  =  SCTDIM  *  RAN(SEEDl)  *  SIGN1 

YRAND  =  SCTDIM  *  RAN(SEED2)  *  SIGN2 

SEED1  =  SEED1  +  .13 

SEED2  =  SEED2  +  .133 

SEED3  =  SEED3  +  .1333 

SEED4  =  SEED4  +  .13333 

TEMP3  =  XRAND/ ASFSCT 

TEMP4  =  YRAND/ ASFSCT 

TEMP 2  =  ASFSCT*SQRT( ( TEMP 3+XI 1-CENTER) **2  +  (TEMP4+XJ1-CENTER)**2 ) 

IF(TEMP2.GT. SCTRAD)  GOTO  1000 

NUMBER  =  NUMBER  +  1 

XLEG2  =  SQRT(XRAND**2  +  YRAND**2) 

HYPO  =  SQRT(XLEG2**2  +  RANGE**2) 

C  =  ( SCAT/ HYPO )*CEXP( (0.0,1 .0)*XK*(RANGE+HYP0) )  +  C 
IF( NUMBER. EQ.NUM)  GOTO  1004 
1000  CONTINUE 

1004  PRS(Il.Jl)  =  C 

1005  CONTINUE 
1010  CONTINUE 

RETURN 
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C*  *  **  *  *  *  *  **  *  *  **  *  *  *  *  **  *  *  *  *  *  *  *  **  *  ***  **  *  *  *  **  **  *  *  *  *  *  *  *  ***  *  Jit  **  **  **  **  **  *  *  *  *  **  *  * 

c 

C  VOLSCT:  THIS  SUBROUTINE  RANDOMLY  PLACES  A  PARTICLE  WITHIN  A  CYLINDER 
C 

C* ************************************************************* ********** 

SUBROUTINE  VOLSCT 
COMPLEX  C,PRS(4i,41) 

PCT  =  FLOAT(NUM) 

SEEDI  =  .012345 
SEED2  =  .543210 
SEED3  =  .975310 
SEED4  =  .013579 
SEEDS  =  .024531 
DO  2010  II  =  10,32 
DO  2005  J1  =  10,32 

TEMPI  =  ASF  *  SQRT( (FLOAT( II )-CENTER)**2+( FLOAT(Jl)- CENTER)** 2) 
IF(TEMPl.GT.TRANR)  GOTO  2005 
NUMBER  =  0 
C  =  (0.0, 0.0) 

PIC(Il.Jl)  =  'X' 

DO  2000  12  =1,999999 

SIGN1  =1.0 

SIGN2  =1.0 

TEMP 5  =  RAN( SEED3) 

TEMP6  =  RAN(SEED4) 

IF( TEMP5 .LT. .5000)  SIGN1  =  -1.0 

IF( TEMP6.LT. .5000)  SIGN2  =  -1.0 

XRAND  =  APERT  *  RAN( SEEDI)  *  SIGN1 

YRAND  =  APERT  *  RAN(SEED2)  *  SIGN2 

SEEDI  =  SEEDI  +  .13 

SEED2  =  SEED2  +  .133 

SEED3  =  SEED3  +  .1333 

SEED4  =  SEED4  +  .13333 

SEED5  =  SEED 5  +  .133333 

TEMP 3  =  XRAND/ ASF 

XII  =  FLOAT(Il) 

XJ1  =  FLOAT(Jl) 

TEMP 4  =  YRAND/ ASF 

TEMP 2  =  ASF* SQRT((TEMP3+XI1 -CENTER) **2+(TEMP4+XJ 1-CENTER) **2) 

IF(TEMP2.GT.TRANR)  GOTO  2000 

NUMBER  =  NUMBER  +  1 

XLEG2  =  SQRT(XRAND**2  +  YRAND**2) 

DEPTH  =  (RAN(SEED5)  *  .1)  +  RANGE 
HYPO  =  SQRT(XLEG2**2  +  DEPTH**2) 

C  =  ( SCAT/ HYPO)* CEXP( (0.0, 1 .0)*XK*(DEPTH+HYPO) )  +  C 
IF( NUMBER. EQ.NUM)  GOTO  2004 
2000  CONTINUE 

2004  PRS(I1,J1)  =  C 

2005  CONTINUE 
2010  CONTINUE 

RETURN 

END 
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Q********** *********** ** *************************************************** 

c 

C  RDNSCT:  THIS  SUBROUTINE  TAKES  THE  COORDINATE  II, J1  FROM  THE  TRANSDUCER 
C  FACE  AND  SUMS  THE  AMPLITUDE  AND  PHASE  DUE  TO  EACH  SCATTERER  WITHIN 
C  THE  BEAM  PROFILE.  ASSUME  A  PLANAR  INCIDENT  WAVE  WITH  A  SPHERICAL 
C  REFLECTED  WAVEFRONT 

C  1)  ASF  -  ARRAY  SCALE  FACTOR  (CM/ //BLOCKS  FOR  TRANSDUCER) 

C  2)  PRS  =  COMPLEX  PRESSURE  ARRAY 

C  3)  XRAND  =  RANDOM  X  COORDINATE  (YRAND  -  Y  COORDINATE) 

C 

£************************************************************************* 

SUBROUTINE  RDNSCT 
COMPLEX  C,PRS(41,41) 

PCT  =  FLOAT(NUM) 

SEED1  =  .012345 
SEED2  -  .543210 
SEED3  =  .975310 
SEED4  =  .013579 
DO  3010  II  =■  l.M 
DO  3005  J1  =  l.M 
C  -  (0.0, 0.0) 

NUMBER  «  0 
XII  =■  FLOAT(Il) 

XJ1  -  FLOAT(Jl) 

TEMPI  =  ASF*SQRT( (XI1-CENTER)**2  +  (XJ 1-CENTER) **2) 

IF( TEMPI. GT.TRANR)  GOTO  3005 
PIC(I1,J1)  -  'X' 

DO  3000  12  =  1,10000 
SIGN1  *  1.0 
SIGN2  ■  1.0 
TEMP5  =*  RAN(  SEED3) 

TEMP 6  =  RAN(SEED4) 

IF( TEMP5.LT. .5000)  SIGNl’-  -1.0 

IF( TEMP6.LT. .5000)  SIGN2  =  -1.0 

XRAND  =  APERT  *  RAN(SEEDl)  *  SIGN1 

YRAND  »  APERT  *  RAN(SEED2)  *  SIGN2 

SEED1  =*  SEED1  +  .13 

SEED2  =■  SEED  2  +  .133 

SEED3  =  SEED 3  +  .1333 

SEED4  -  SEED4  +  .13333 

TEMP3  =■  XRAND/ ASF 

TEMP 4  =  YRAND/ ASF 

TEMP2  =■  ASF*SQRT( ( TEMP3+XI l-CENTER)**2+( TEMP4+XJ 1-CENTER) **2 ) 

IF(TEMP2. GT.TRANR)  GOTO  3000 

NUMBER  =*  NUMBER  +  1 

XLEG2  =  SQRT(XRAND**2  +  YRAND**2) 

HYPO  =  SQRT(XLEG2**2  +  RANGE**2) 

C  -  ( SCAT/ HYPO ) * CEXP( (0.0, 1.0) *XK* ( RANGE+HYPO ) )  +  C 
IF( NUMBER. EQ.NUM)  GOTO  3004 
3000  CONTINUE 

3004  PRS(Il.Jl)  =  C 

3005  CONTINUE 
3010  CONTINUE 


n  n  n  n  n  no 
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Q* ********************************************************************* 

FILINC:  INCREMENT  DATA  FILE  NAMES  &  CALL  INFO  AND  SAVDAT  SUBROUTINES 


SUBROUTINE  FILINC 

DOUBLE  PRECISION  FNAMEI , FNAME2 , FNAME3 , FNAME4 

FLAG  OPERATIONS: 

0  =  PERFORM  ROUTINE 
1  -  SKIP 

IF(IFLAG.EQ.l.AND.IFLAGl.EQ.l)  GOTO  4005 

FNAMEI  -  'PLT000.DAT' 

FNAME2  =  'AMP000.DAT' 

FNAME3  =  ' INF000.DAT' 

FNAME4  =  'PHZ000.DAT' 

C  ALLOW  INCREMENT  OF  FILE  NAMES  DEPENDING  UPON  THE  CHANGING  VARIABLE 

IF(ANS.EQ.'A')  ICOUNT  -  IFIX( APERT  *  100.00) 

IF(ANS.EQ.'R')  ICOUNT  *  IF IX( RANGE  *  10.00) 

IF(ANS.EQ.'F')  ICOUNT  *  IFIX(FREQ  *  10.00) 

IF(ANS.EQ.'N')  ICOUNT  *  NUM 

DO  4000  I  =■  1,  ICOUNT 
CALL  INCFIL( FNAMEI) 

CALL  INCFIL(FNAME2) 

CALL  INCFIL(FNAME3) 

4000  CALL  INCFIL( FNAME4 ) 

IF(IFLAG.EQ.O)  CALL  PLTDAT( FNAMEI ) 

IF(IFLAGl.EQ.O)  CALL  PRTDAT( FNAME2 , FNAME4 ) 

IF(IFLAGl.EQ.O)  CALL  INFO ( FNAME 1 , FNAME3 ) 

4005  RETURN 
END 


o  o  o 
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c********************************************* ************************** 
INCVAR:  INCREMENT  COUNTER  VARIABLES  AND  SET  CONTROL  FLAG  2 
C*********************************************************************** 


SUBROUTINE  INCVAR 
IFLAG2  -  0 

IF(ANS.NE.'N')  GOTO  5000 
NUM  =  NUM  +  IFIX(XINC) 

ISTOP  -  IFIX(STOP) 

IF(NUM.GT.ISTOP)  IFLAG2  =*  I 
GOTO  5015 

5000  IF(ANS.EQ.'A' .OR.ANS.EQ.'F')  GOTO  5005 
RANGE  *  RANGE  +  XINC 
IF( RANGE. GT. STOP)  IFLAG2  *  1 
GOTO  5015 

5005  IF(ANS.EQ.'F')  GOTO  5010 
APERT  -  APERT  +  XINC 
IF  (  APERT .  GT .  STOP )  IFLAG2  **  1 
GOTO  5015 

5010  FREQ  -  FREQ  +  XINC 

IF( FREQ. GT.  STOP)  IFLAG2  =*  1 

5015  RETURN 
END 


o  o 
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Q* **************************** **************************************** 

c 

C  LINSCT:  THIS  SUBROUTINE  TAKES  EACH  COORDINATE  ON  THE  TRANSDUCER 
C  AND  SUMS  THE  AMPLITUDE  AND  PHASE  DISTRIBUTION  DUE  TO  EACH  SCATTER. 
THE  SCATTERERS  ARE  ARRANGED  IN  A  LINEAR  (ONE-DIMENSIONAL)  ARRAY. 

c*  ******************************************************************** 


SUBROUTINE  LINSCT 
COMPLEX  C, PRS(4 1 ,41) 

XEDGE  =  IFIX( CENTER)  -  INT(DIM/2.0) 

YEDGE  =  CENTER 
PCT  =  FLOAT(NUM) 

DO  6020  I  =1 ,M 
DO  6015  J  -1,M 
XI  =  FLOAT(I) 

XJ  =  FLOAT(J) 

XTEMP  =  ASF*SQRT( (X1-CENTER)**2  +  (XJ-CENTER)**2) 
IF(XTEMP.GT.TRANR)  GOTO  6015 
C  =  (0.0, 0.0) 

PIC(I,J)  =  'X' 

IF(NUM.NE.l)  GOTO  6000 

XLEG2  =  SQRT( ( ( ( ASF ) *(XI-XEDGE )- ( TRANR) ) **2 )  + 

1  ( ( (XJ-YEDGE)*(ASF) )**2) ) 

HYPO  *  SQRT(XLEG2**2  +  RANGE**2) 

C  =( SCAT/ HYPO )*CEXP((0. 0,1.0)  *XK* ( RANGE+HYPO ) ) 

GOTO  6011 

6000  DO  6010  12  =1,NUM 
XI2  -  FLOAT( 12) 

XLEG2  =SQRT(((ASF)*(XI-XEDGE)-(XI2-1.0)*(PSF))**2  + 

1  ( (XJ-YEDGE)*( ASF) )**2) 

HYPO  «■  SQRT(XLEG2**2  +  RANGE**2) 

C  =  ( SCAT/ HYPO )*CEXP( (0.0, 1.0)  *  XK  *  (RANGE+HYPO))  +  C 

6010  CONTINUE 

6011  PRS( I, J )  =  C 
6015  CONTINUE 
6020  CONTINUE 

RETURN 

END 


, 
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C* ******************************************************************** 

C 

C  PARCNT :  LOADS  PARTICLE  COUNT  ARRAY  (RPCT)  USED  IN  THE  RECTANGULAR 
C  SCATTERERING  ARRANGEMENT  IN  THE  NUMBER  COUNT  ARRAY  (RNUM)  FOR  USE  IN 
C  SAVDAT  SUBROUTINE. 

C 

£* ************************************************************* ******* 


SUBROUTINE  PARCNT 


C  SCTANS : 

C  L,R,V,A  =  EXIT 
CO  =  PERFORM  TRANSFER 

T  *  SCTANS 

IF^.EQ.'L'.OR.T.EQ.'R'.OR.T.EQ.'V'.OR.T.EQ.'A')  GOTO  7010 

DO  7000  I  =  1 , NCOUNT 
RNUM(I)  =  RPCT(I) 

7000  CONTINUE 


TYPE  7005 

7005  FORMAT ( ' ODATA  TRANSFER  FROM  RPCT  ARRAY  TO  RNUM  COMPLETED ! ' ) 

7010  RETURN 
END 

c************************************************************************* 

C 

C  DECIDE:  THIS  SUBROUTINE  CALLS  THE  CORRECT  SCATTERERING  ROUTINE  BASED 
C  ON  USER  INPUT. 

C 

C* ******* ***************************************** ************************ 

SUBROUTINE  DECIDE 


C  SCATTERER  ROUTINE:  USER  DEFINED 
C  L  =  LINEAR  ARRAY 

C  0  =  RECTANGULAR  ARRAY 

C  R  =  RANDOM  ARRAY 

C  V  =  RANDOM  VOLUME 

C  A  *  RANDOM  ARRAY  CONSTANT  SCATTERING  APERTURE 


IF( SCTANS. EQ.'L') 
IF( SCTANS. EQ. ' O' ) 
IF( SCTANS. EQ.'R' ) 
IF( SCTANS. EQ.'V') 
IF( SCTANS. EQ.' A') 
RETURN 
END 


CALL  LINSCT 
CALL  RCTSCT 
CALL  RDNSCT 
CALL  VOLSCT 
CALL  RDNAPT 


o  o  o  o 
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C* ***************************************** ************ *************** 


RCTSCT :  BACKS CATTER  ROUTINE  FEATURING  A  RECTANGULAR  ARANGEMENT 
OF  PARTICLES  WITHIN  A  TRANSDUCER  BEAM  PROFILE 

C* ********************** ********************************************** 


SUBROUTINE  RCTSCT 
COMPLEX  C, PRS(4 I ,41) 

XEDGE  =  10.0 

YEDGE  =  10.0 

TEMP 5  =  FL0AT(NUM)/2.0 

TEMP 2  =*  TEMP 5  -  INT(TEMP5) 

IF( TEMP 2 . NE . 0 . 00 )  GOTO  3000 
PCNTR  =  TEMP5  +  .5000 
GOTO  3005 

3000  PCNTR  =  INT(TEMP5)  +1.0 

3005  DO  3025  II  =  10,32 

DO  3020  J1  =  10,32 

XII  =  FLOAT(Il) 

XJ1  =  FLOAT(Jl) 

TEMPI  =  ASF*SQRT( (XI 1-CENTER) **2  +  (XJ1-CENTER)**2) 

IF ( TEMP l.GT. TRANR)  GOTO  3020 
C  =  (0.0, 0.0) 

PIC(I1,J1)  =  'X' 

DO  3015  12  =  1 , NUM 
DO  3010  J2  =  l.NUM 

TEMP3=PSF*SQRT((FLOAT(I2)-PCNTR)**2+(FLOAT(J2)-PCNTR)**2) 

IF ( TEMP3 . GT . TRANR)  GOTO  3010 

XLEG2=SQRT( ( ( ASF)*(XI1 -XEDGE )- ( FLOAT( 12 )- 1 . 0 )* ( PSF) ) **2  + 

1  ((XJ1-YEDGE)*(ASF)-(FL0AT(J2)-1.0)*(PSF))**2) 

HYPO  =  SQRT(XLEG2**2  +  RANGE* *2) 

C  »( SCAT/ HYPO )*CEXP( (0.0, 1 .0)*XK*(RANGE+HYPO) )  +  C 
IF(XJ1. NE. CENTER. OR.XI1.NE. CENTER)  GOTO  3010 
IF( 12 -NE. IF IX( PCNTR) .OR. J2. NE. IF IX( PCNTR) )  GOTO  3010 
PCT  »  0.0 

DO  3007  13  =  l.NUM 
DO  3006  J3  =  1 ,NUM 

TEMP4=PSF*SQRT( (FLOAT( I3)-PCNTR)**2  +  (FLOAT(J3)-PCNTR)**2) 
IF(TEMP4.LE. TRANR)  PCT  =  PCT  +  1 

3006  CONTINUE 

3007  CONTINUE 
3010  CONTINUE 
3015  CONTINUE 

PRS(Il.Jl)  =  C 
3020  CONTINUE 
3025  CONTINUE 
RETURN 
END 


o  o  o 


C**************** ******************************************** *********** 

c 

C  WARN:  CHECK  PARAMETERS  TO  SEE  IF  DATA  ASSUMPTIONS  HOLD 
C 

Q*********************************************************************** 


SUBROUTINE  WARN 

4003  XKA  =  XK  *  PARRAD  *  1.0E-04 
IF(XKA.GT.l.O)  TYPE  4005 

4005  FORMAT( ' ***KA  >  1.0  CHECK  INITIAL  CONDITIONS***') 
IF( RANGE. LT.XNEAR)  TYPE  4010 

4010  FORMAT('0***SCATTERERING  ARRAY  WITHIN  NEARFIELD***') 
RETURN 
END 


C***  *  ****************************************************************** 


VOLDAT:  WRITE  AVERAGE  RECEIVED  PRESSURE  PER  VOLUME  TO  DISK 

Q* * ****** **★★★*★**★ ************ ********************************* ******* 


SUBROUTINE  VOLDAT 
DOUBLE  PRECISION  FNAME 

IF(SCTANS.NE.'V')  GOTO  5015 
FNAME  *  'VD000.DAT' 

ICOUNT  =  IF IX( APERT  *  100.00) 

DO  5000  I  =1, ICOUNT 
5000  CALL  INCFILC FNAME) 

CALL  OPEN(21, 'FILE' , FNAME) 

DO  5005  I  =1 , NCOUNT 
WRITE(21 , 5010)  RNUM( I ) , RVOL( I ) 
5005  CONTINUE 

5010  F0RMAT(1X,F11.6,1X,F11.6) 

CALL  CL0SE(21) 

5015  RETURN 
END 


o  u  o 
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C*  *  *  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  **  *  *  *  *  *  **  *  **  *  *  **  *  *  *  *  *  **  *  *  *  *  *  *  *  **  *  *  **  **  * 

ASKFEW:  THIS  SUBROUTINE  IS  USED  WHEN  A  MACRO  FILE  CONTROLS  EXECUTION 

C* ************************************************* ******************** 

SUBROUTINE  ASKFEW 

ACCEPT  8000, NUM 
8000  FORMAT(I) 

8005  FORMAT(F) 

XINC  =  1.0 
STOP  =  100.0 
APERT  =  .25 
FREQ  =«  5.0 
RANGE  =  15.0 
IFLAG1  =  0 
ANS  =  'N' 

IF LAG  =  1 
SCTANS  =  'V' 

RETURN 

END 


C******************************************************************** 

c 

c  APPSTR:  APPENDS  ONE  CHARACTER  STRING  TO  ANOTHER 
C  1)  STRING  -  STRING  TO  BE  EXTENDED 

C  2)  NCH  -  NUMBER  OF  CHARACTERS  IN  STRING  UPON  ENTRY 

C  3)  ADDTN  -  STRING  TO  BE  APPENDED 

C  4)  NADD  -  NUMBER  OF  CHARACTERS  TO  BE  APPENDED 

C  5)  167  -  CHARACTER  TYPE  :  6  -  6-BIT  (FROM  ASCPK) 

C  7  =  7-BIT  (STD.  ASCII) 

C 

Q******************************************************************** 

SUBROUTINE  APPSTR( STRING, NCH, ADDTN, NADD, 167 ) 

DIMENSION  STRING(l) 

NBWD  =■  (36/167)  *  167 

NBITS  =  NCH  *  167 

NWD  -  NBITS/ NBWD 

NBITR  =*  NBITS  -  NWD  *  NBWD  -  1 

NWD  =*  NWD  +  1 

CALL  BPOINT( 167 , STRING( NWD) , NBITR, BPS ) 

CALL  BP0INT( 167 , ADDTN, -1 , BPA) 

DO  9000  I  =1 ,NADD 
CALL  ILDB( ICH, BPA) 

CALL  IDPB( ICH, BPS) 

9000  CONTINUE 


RETURN 

END 
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C********************************************************************** 

STAT:  CALCULATE  USING  AMP(I,J)  AND  THETA(I, J)  ARRAYS  STATISTICAL 
PARAMETERS  -  MEAN  AND  VARIANCE. 

C********************************************************************** 

SUBROUTINE  STAT 

K  -  0 

TEMPA  «  0.0 
TEMPB  =*0.0 
DO  1005  I  -  1,M 

DO  1000  J  =  1,M 
IF(PIC(I,J).EQ.'-')  GOTO  1000 
TEMPA  =  TEMPA  +  AMP(I,J) 

TEMPB  -  TEMPB  +  THETA(I.J) 

K  =*  K  +  1 

1000  CONTINUE 
1005  CONTINUE 

AMEAN  -  TEMPA/ FLOAT(K) 

PMEAN  =*  TEMPB/ FLOAT(K) 

K  -  0 

TEMPA  -  0.0 
TEMPB  -  0.0 
DO  1015  I  -  1,M 

DO  1010  J  **  1  ,M 
IF(PIC(I,J).EQ.'-')  GOTO  1010 
TEMPA  =  TEMPA  +  (AMP(I.J)  -  AMEAN)** 2 
TEMPB  =  TEMPB  +  (THETA( I, J)  -  PMEAN)**2 
K  =  K  +  1 

CONTINUE 
CONTINUE 

AVAR  -  TEMPA/ (FLOAT(K)  -  1.0) 

C  PVAR  *  TEMPB/ (FLOAT(K)  -  1.0) 

TYPE  1020, AMEAN, PMEAN 

1020  FORMAT( ' OMEAN  AMPLITUDE  -  ',F11.6,'  MEAN  PHASE  -  ',F11.6) 

TYPE  102 5, AVAR 

1025  FORMAT('OAMP  VARIANCE  =',F11.6) 


C 

1010 

1015 


RETURN 

END 
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C***************************** ******************************************* 


NAME:  CONSTRUCTS  OUTPUT  DATA  FILE  NAME  USING  USER  SPECIFIED  PARAMETERS 
FIRST  CHARACTER  REFERS  TO  THE  SCATTERING  ARRANGEMENT  AND  THE 
SECOND  INDICATES  WHICH  PARAMETER  IS  BEING  VARIED 


C***************  ********************************************************* 


SUBROUTINE  NAME( FNAME5 , FNAME6 ) 

DOUBLE  PRECISION  FI , FNAME5 , FNAME6, S2 

A2  -  ANS 
S2  -  SCTANS 
FI  -  '0000.DAT' 

FNAME6  *  'TAB000.DAT' 

INAME  *  IF IX( APERT  *  100.00) 

CALL  APPSTR(S2,1,A2, 1,7) 

CALL  APPSTR(S2,2,F1,4,7) 

FNAME5  =  S2 

DO  2000  I  =  1, INAME 
CALL  INCFIL(FNAME5) 

CALL  INCFIL(FNAME6) 

2000  CONTINUE 

RETURN 

END 
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PLT3D 

The  graphics  program,  PLT3D,  plots  the  three-dimensional  ampli¬ 
tude  and  phase  distributions  across  a  transducer  aperture.  This  pro¬ 
gram  can  be  executed  at  installations  supporting  the  TOPS- 10  operating 
system  and  Tektronix  AG-II  graphics  software.  PLT3D  operates  in  an 
interactive  mode  allowing  flexible  control  of  graphics  parameters  and 
subsequent  plot  generation.  The  main  graphics  subroutine,  VISDO  (Visu¬ 
alize  Double  Surfaces),  requires  data  to  be  expressed  in  cartesian 
coordinates.  VISDO  plots  an  array  of  M  x  N  elements,  each  array  ele¬ 
ment  having  a  distinct  value.  The  user  specifies  the  shear  or  rotation, 
relative  elevation  and  span  of  the  plot.  The  user  may  also  specify 
the  plot  not  to  be  cross-hatched,  by  including  the  subroutine  TSPLT, 
meaning  only  lines  in  the  left-to- right  direction  are  drawn  rather 
than  left-to-right  and  front-to-back.  The  algorithm  used  by  TSPLT 
accounts  for  hidden  lines  thus  maintaining  data  integrity. 


CPLOT 

The  FORTRAN  program  CPLOT  constructs  contour  plots  using  data 
generated  by  PHASE.  Once  the  amplitude  and  phase  data  have  been  read 
into  program  arrays,  the  number  of  contour  levels  must  be  specified. 
This  parameter  is  used  to  determine  how  many  contours  of  equivalue 
are  drawn  between  the  maximum  and  minimum  data  points.  The  contour 
subroutine  PLTKP  then  drawn  the  final  plot.  This  plot  only  shows  the 
regularity  of  the  data  in  one  plane  with  no  evidence  of  relative 
height  between  contours. 


n  n 
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(2****************  it  it  irk  it  it  irk  it  irk  ******  it  irk  it*  it  ******  it  irk  ******  it  it  Irk  it  ******* 


PHASE-CANCELLATION  STUDIES:  3D  PLOTTING  PROGRAM 
C 

C  THIS  PROGRAM  USES  DATA  FILES  CREATED  BY  PHASE. F4  AND  CONSTRUCTS 
C  3-D  PLOTS.  THE  PROGRAM  REQUESTS  U  AND  V  SHEAR  DIRECTIONS 
C  EAST  OR  WEST  ORIENTATION  AND  HAS  FLEXIBLE  CONTROL  OVER  SCALE 
C  AND  COMPRESSION  FACTORS. 

C 

C  RUN  COMMAND :  PLT3D . F4 , NUMBER . REL , @  SYS : GRA3D , 0SYS : GRALIB 
C 

c* ************************ ************************************* **** 

DOUELE  PRECISION  FNAME1 
D IKC:: SIC*!  AMP ( 4 1 , 4 1 )  , THETA( 4 1 , 4 1 ) 

EXTERNAL  PLTCA 

TYPE  1001 

1001  F0RMAT('0  ') 

TYPE  1002 

TYPE  1001 
1000  TYPE  1005 

1002  FORMAT(20X,' BIOENGINEERING:  3-D  GRAPHICS  PROGRAM') 

1005  FORMATC'O ENTER  DATA  FILE  NAME:  ',$) 

ACCEPT  1010 ,FNAME1 
1010  FORMAT(AIO) 

C  READ  IN  DATA  FILE  CREATED  BY  PHASE. F4 
CALL  0PEN(21 , ' FILE' ,FNAME1 ) 

READ( 21,1025)  M, AMPMIN , AMPMAX , THEMIN, THEMAX , SCTANS 

1025  FORMAT( 13, IX, FI  1.5, IX, FI 1.5, IX, FI  1.5, IX, FI  1.5, IX, Al) 

READ( 21,1026)  HUM , TRANR, RANGE , FREQ , APERT , PCT 

1026  FORMAT( 15, IX, FI  1.5, IX, FI 1.5, IX, FI  1.5, IX, FI  1.5, IX, FI  1.5) 

READ(  21 , 1027 )  (  (  AMP(  I , J)  ,  J=1  ,M)  ,  I=*l  ,M) 

READ(21 , 1027) ( (THETA( I,J) ,J=1 ,M) ,1=1 ,M) 

1027  F0RMAT(41(F1 1 .6 , IX) ) 

CALL  CLOSE( 21) 

CALL  BELL 

C  MULTIPLY  AMPLITUDE  AND  PHASE  DATA  BY  XSIGN  TO  ACCOUNT  FOR  PLOTTING 
C  BY  SUBROUTINE  VISDO 

XSIGN  =  -1.0 

DO  1039  I  =1  ,M 

DO  1038  J  =1 ,M 

AMP( I,J)  =  XSIGN  *  AMP( I, J) 

1038  THETA(  I ,  J)  =  XSICN  *  TllETA(I.J) 

1039  CONTINUE 
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C  PROMPT  USER  FOR  PLOTTING  NAME  AND  SCALE  FACTORS  TO  BE  USED  BY  VISDO 
1040  TYPE  1045 

1045  FORMAT( 'OENTER  PLOT  NAME  (5  CHARACTERS):  ' ,$) 

ACCEPT  1050 ,PNAME 
1050  FORMAT(A5) 

TYPE  1055 

1055  FORMAT( 'OEAST/WEST  VIEW  (EAST  =  -l.WEST  =  +1):  ' ,$) 

ACCEPT  1060.LVIEW 
1060  FORMAT(I) 

TYPE  1065 

1065  FORMAT( ' 0E1ITER  U  SHEAR  (REAL  0.0  TO  1.0):  '  ,$) 

ACCEPT  1070.USHEAR 
1070  FORMAT(F) 

TYPE  1075 

1075  FORMAT( ' OENTER  V  SHEAR  (REAL  0.0  -  1.0):  '  ,$) 

ACCEPT  1070, VSHEAR 
TYPE  1080 

1080  FORMAT( 'OENTER  COMPRESSION  FACTOR:  ',$) 

ACCEPT  1070, COMP 
TYPE  1085 

1035  FORMAT( 'OENTER  SPAN  COMPRESSION  FACTOR:  ',$) 

ACCEPT  1070, SFACTR 
ASP  AN  1  =■  1.0  *  COMP 
ASPAN2  =  -ASPAN1 
ASPAN1  =  ASPAN1/ SFACTR 
TSPAN1  =  (THEMAX  -  THEMIN)*COMP 
TSPAN2  *  -TSPAN1 
TSPAN2  =  TSPAN2/ SFACTR 
1141  TYPE  1145 

1145  FORMAT('ODO  YOU  WISH  TO  PLOT  AMPLITUDE  OR  PHASE  (A/P):  ',$) 
ACCEPT  1150, ANS 

IF( ANS.NE. ' A' .AND.ANS.NE.'P')  GOTO  1141 

C  SET  UP  GRAPHIC  LIBRARY  AND  CALL  PLOTTING  ROUTINES 

CALL  PLTOO 
CALL  BELL 
CALL  FRAME 

CALL  PRT( ANS , PCT , TRANR , RANGE , FREQ , SCTANS) 

CALL  PLTL.U  PNA11E) 

1150  FORMAT(Al) 

IF(ANS.EQ.'P')  COTO  1155 

CALL  VISDO( ASPAN1 , AMP ,AMP , ASPAN2 ,M,M,M,M,USHEAR, VSHEAR, 

1  LVIEW.-l ,PLTCA) 

GOTO  1160 

1155  CALL  VISD0(TSPAN1, THETA, THETA, TSPAN2 ,M,M,M,M,USHEAR, VSHEAR, 

1  LVIEW.-l , PLTCA) 

1160  CONTINUE 
CALL  BELL 
STOP 
END 
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C* ************************************** *************************** 

C 

C  PRT:  PRINT  NUMERICAL  INFORMATION  ON  GRAPH 
C 

(^****************************************************************** 


SUBROUTINE  PRT( ANS , PCT , TRANR , RANGE , FREQ , SCTANS) 

DIMENSION  A(6),B(5),C(5),D(2),EC2),F(3),G(4),HC3),IC3) 
DATA(A(K),K= 1,6) /'NORMALIZED  RECEIVED  AMPLITUDE'/ 
DATA(3(K),K=1,5)/' NORMALIZED  RECEIVED  PHASE'/ 

DATA(D(K) ,K=1 ,2) /' APERTURE :' / 

DATA( E( K) , K* 1 , 2 ) / ' DISTANCE : ' / 

DATA( F( !C), K=  1,3) /'RANDOM  ARRAY'/ 

DATA( G(K) ,K=1 ,4)/' RECTANCULAR  ARRAY ' / 

DATA(H(K)  ,K-1 ,3)/ 'LINEAR  ARRAY'/ 

DATA( I( K) ,K-1 , 3)/' RANDOM  VOLUME'/ 

APERT  -TRANR  *  2.0 
CALL  CHRSET(21 ,33) 

IF(ANS.EQ.'P')  GOTO  2000 
CALL  JUSTFX(  30 , A,0  ,LEN,  IOFS) 

CALL  MOVABS( 5 1 2+IOFS , 725) 

CALL  ANCHOS(A.LEN) 

GOTO  2005 

2000  CALL  JUSTFX( 25 , B,0 ,LEN, IOFS) 

CALL  MOV ABS( 51 2+IOFS, 725) 

CALL  ANCHOS(B.LEN) 

2005  CALL  CHRSET( 14,22) 

CALL  MOVABSC 100,655) 

STR  * ' NUM : ' 

NUM  -  IFIX(PCT) 

CALL  ANCHOS( STR, 4) 

CALL  NUMBER( NUM , ' 15 ' ) 

CALL  MOVABS( 100,630) 

STRING  3 ' CM' 

CALL  JUSTFX( 10, D, -l ,LEN,IOFS) 

CALL  ANCHOS(D,LEN) 

CALL  NUMBER( APERT , ' F6 . 2 ' ) 

CALL  ANCHOSC STRING, 2) 

CALL  MOV ABS( 100,605) 

CALL  ANCHOS(E,9) 

CALL  NUMBER( RANGE, 'F7 .2' ) 

CALL  ANC110S( STRING, 2) 

CALL  MOVABS( 100,580) 

STR  -'FREQ:' 

CALL  ANCHOSC STR, 5) 

CALL  NUMBER( FREQ , ' F  5 . 2 ' ) 

STR  -'MHZ' 

CALL  ANCHOSC STR, 3) 

CALL  MOVABSC 100,555) 

IF C SCTANS. EQ.'R' .OR. SCTANS .EQ. ' L' .OR. SCTANS. EQ. 'V' )  GOTO  2010 
CALL  JUSTFXC 1 7 ,C,-1 , LEN, IOFS) 

CALL  ANCHOSC G, LEN) 

GOTO  2030 


94 


2010  IF ( SCTANS . EQ . ' R' . OR. SCTANS •  EQ . ' V ' )  GOTO  2020 

CALL  JUSTFX(  12 ,11,-1  , LEN, IOFS) 

CALL  ANCHOS(H.LEN) 

GOTO  2030 

2020  IF( SCTANS. EQ.'R')  COTO  2025 

CALL  JUSTFX( 13, 1,-1 ,  LEN, IOFS) 

CALL  ANCHOS(I.LEN) 

GOTO  2030 

2025  CALL  JUSTFX( 12.F.-1 , LEN, IOFS) 

CALL  ANCHOS(F.LEN) 

2030  RETURN 
END 


Qic  ************** *** *********** ******************** ******* ********* ****** 

c 

C  FRAME:  DRAW  A  FRAME  AROUND  THE  PLOT.  SET  UP  AN  8.5"  X  11"  FRAME 
C  INNER  FRAME  OF  1"  INSIDE  OF  IT  AND  DEFINE  THE  ORIGIN  PAGE  CENTER 
C 

c*********************************************************************** 


SUBROUTINE  FRAME 
CALL  PLOT(0 .0 ,0 .0 , 3) 
CALL  PL0T(0 .0 , 11.0,2) 
CALL  PLQT(3 .5 , 11.0,1) 
CALL  PLOT(8 .5 ,0 .0 ,  1 ) 
CALL  PL0T(0 .0 ,0 .0 , 1) 
CALL  PLOT( 4.25,5.50,-3) 
RETURN 
END 
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C* *********** ******** ****** ************ ******************************* 

C 

C  TSPLT:  PLOT  3-D  REPRESENTATION  ON  DATA  WITHOUT  CROSSHATCHING 
C  U(K)  -  ARRAY  FOR  DISTANCE  ALONG  X  COORDINATE 

C  V(K)  -  ARRAY  CONTAINING  AMPLITUDE  AT  DISTANCE  X 

C  IX(J,I)  -  CONVERTS  2  DIMENSIONAL  ARRAY  INTO  1 
C 

c* ********************************************** ********************** 


SUBROUTINE  TSPLT( ASP AN 1 , AMP , ASPAN2 , N , M , USHEAR , VSHEAR , LV ) 
EXTERNAL  PLTCA 

DIMENSION  AMP( 1 ) ,U( 501 ) ,V(501) 

C  FUNCTION  STATEMENTS  &  CALL  NULL  INITIAL  HORIZON  FOR  HIDDEN  LINE 

IX(J,I)  =  (I  -  1)  *  M  +  J 

AM( J , I)  =  ZS  *  ( AMP( IX( J , I) )  -  ASPAN1 ) 

CALL  VISNH 

C  SET  UP  SCALE  FACTORS  USING  SHEAR  VALUES 
MK  =  501 

ZS  =  (1.0  -  VSHEAR) /(ASPAN2  -  ASPAN1) 

DELTAU  =  USHEAR/ FLOAT(M  -  1) 

DELTAV  =  VSHEAR/ FLO AT(M  -  1) 

DELTAS  =  (1.0  -  USHEAR) /FLOAT(M  -  1) 

VE  =  0.0 

C  SET  UP  DATA  FOR  VISUALIZE  HORIZON  ROUTINE  VISHO 

DO  1010  I  =  1  ,M 
K  =  0 

EU  =  USHEAR  -  (DELTAU-  *  FLOAT(I)) 

DO  1000  J  =  1,M 
K  =  MIN0(K+1 ,MK) 

U(K)  =  EU 

V(K)  =  VE  +  AM( J , I) 

EU  =  EU  +  DELTAS 
1000  CONTINUE 

CALL  VISHO(U,V,K,l , PLTCA) 

CALL  VISHO(U,V,K,-l , PLTCA) 

VE  =  VE  +  DELTAV 
1010  CONTINUE 
RETURN 
END 
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£********************************************************************** 

c 

C  PHASE-CANCELLATION  STUDIES:  CONTOUR  PLOTTINC  PROCRAM 
C 

C  THIS  PROCRAM  USES  DATA  CREATED  FROM  PHASE. F4  AND  PLOTS  USER 
C  SPECIFIABLE  CONTOUR  INTERVALS.  USES  3DPLOT.F4  AS  THE  MAIN  ROUTINE. 

C 

C  LOAD  CPLOT.  F4  ,@SYS  :CRA3D 
C 

Q* ********************************************************************** 

DOUBLE  PRECISION  FNAME1 
DIMENSION  AMP(41 ,41) ,THETA(4I ,41) 

EXTERNAL  PLTCA 

TYPE  1001 

1000  TYPE  1005 

1001  FORMAT( 20X , ' BIOENGINEERING:  CONTOUR  GRAPHICS  PROGRAM') 

1005  FORMAT( 'OPLEASE  ENTER  PLOTTING  DATA  FILE:  ',$) 

ACCEPT  1010.FNAME1 
1010  FORMAT(AIO) 

C  READ  IN  DATA  FILE  AS  SPECIFIED  BY  USER 
CALL  OPEN( 21,' FILE' , FNAME1 ) 

READ( 21,1015)  M , AMPMIN , AMPMAX , THEMIN , THEMAX , SCTANS 
1015  FORMAT( 13, IX, FI 1. 5, IX, FI 1.5, IX, FI 1.5, IX, FI 1.5, IX, Al) 

READ( 21,1020)  NUM , TRAKR , RANGE , FREQ , APERT , PCT 
1020  FORMAT( I5,1X,F11.5,1X,F11.5,1X,F11.5,1X,F11.5,1X,F11.5) 
READ(21,1025)((AMP(I,J),J*l,M),I=l,M) 

READ(21 , 1025 )( (THETA( I,J) ,J=1 ,M) ,1-1 ,M) 

1025  F0RMAT(41 (FI  1 .6, IX) ) 

CALL  CL0SE(21) 

C  SET  UP  PLOTTING  PARAMETERS 

1050  CINT1  =  AMPMIN 
CINT2  =  AMPMAX 
NUMX  -  M 
NUMY  -  M 
SUBX  -5.0 
SUBY  -  5.0 
TYPE  1055 

1055  FORMAT( 'OENTER  PLOT  NAME  (MAX  5  CHARACTERS):  ',$) 

ACCEPT  1060 , FILNAM 
1060  F0RMAT(A5) 

TYPE  1065 

1065  F0RMAT( 'OHOW  MANY  CONTOUR  LEVELS  ( INTECER) :  ' ,$) 

ACCEPT  1070 .NLEVELS 
1070  FORMAT(I) 

C  CALL  GRAPHICS  LIBRARY  ROUTINES  AND  CALL  CONTOUR  ROUTINE  PLTKP 

CALL  PLTOO 
CALL  PLTFR 


CALL  PLTLA(FILNAM) 

CALL  PLTKP( CINT1 , AMP , CINT2 .NLEVEL, SUBX.M , SUBY ,M, PLTCA) 
CALL  BELL 
TYPE  1035 

1035  FORMAT('OREPEAT(Y/N):  ' ,$) 

ACCEPT  1030.ANS 
IF(ANS.EQ.'Y')  GOTO  1050 
STOP 
END 
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BATGEN 

BATGEN  is  a  FORTRAN  program  written  to  control  the  batch  stream 
of  the  PDP-10  at  the  Hybrid  Computer  Laboratory.  The  program  reads 
a  parameter  data  file  after  each  successful  execution  of  PHASE.  BATGEN 
then  writes  a  new  batch  control  file  and  resubmits  PHASE  to  the  batch 
controller  for  execution.  The  parameter  data  file,  named  PARAM.DAT, 
contains  input  values  of  aperture,  frequency,  range,  and  the  number 
of  scatterers.  This  file  may  also  contain  start,  stop,  and  increment 
values  if  the  user  is  changing  one  of  the  above  parameters.  The  user 
must  ensure  that  disk  memory  allocation  is  not  exceeded  when  using 
the  batch  controller  or  program  execution  will  fail. 


I 
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C* ********************************************************************** 

C 

C  BATCH  CONTROL  FILE  GENERATOR  -  BATGEN.F4 
C 

C  PROGRAM  DESCRIPTION:  PARAM.DAT  CONTAINS  THE  8  INPUTS  TO  PHASE. F4 
C  WHEN  USING  THE  ASKFEW  SUBROUTINE.  THE  PROGRAM  READS  IN  THE  DATA 
C  FILE  AND  REWRITES  IT  TO  A  CONTROL  FILE  NAMED  RDN.CTL  FOR  BATCH  USE. 

C 

c******* ************************** ************************************** 
REAL  PARAM( 10,8) 

CALL  OPEN( 21, 'FILE', 'PARAM.DAT') 

ILINE  =  1 

1 000  READ( 2 1 , 1 005 , END= 1 0 1 0 ) ( PARAM( ILINE , J ) , J=1 , 8 ) 

1005  F0RMAT(8F) 

ILINE  =  ILINE  +  1 
GOTO  1000 

1010  CLOSE(UNIT=21 ,DISPOSE='DELETE' ) 

C  NOW  THAT  THE  DATA  FILE  HAS  BEEN  LOADED  INTO  PARAM  ARRAY,  CHECK 
C  LAST  OPERATION,  REWRITE  DATA  BACK  TO  DISK  OMITTING  FIRST  LINE 


NLINE  =  ILINE  -  1 
IF(NLINE.EQ.l)  GOTO  1025 

CALL  OPEN(21, 'FILE' , 'PARAM. DAT') 

DO  1020  ILINE  =2, NLINE 

WRITE(21 , I015)(PARAM( ILINE, J) ,J*1 ,8) 

1015  FORMAT(3F10.5) 

1020  CONTINUE 

CALL  CLOSE(21) 

C  WRITE  BATCH  CONTROL  FILE,  CHANGE  LOCATION  1  &  5  IN  PARAM  ARRAY  TO 
C  INTEGER  VARIABLES  FOR  INPUT  TO  ASKFEW  SUBROUTINE 

1025  CALL  OPEN( 2 1,' FILE' ,' RDN.CTL') 

I  PARI  =>  PARAM(l.l) 

IPAR5  =  PARAM( 1,5) 

WRITE( 21,1030)  IPAR1 , PARAM( 1,2), PARAM( 1,3), PARAM( 1,4), 

1  IPAR5 , PARAM( 1,6), PARAM( 1,7), PARAM( 1,8) 

1030  FORMAT( ' . SET  TIME  0',/, 

1  '.SET  HPQ  0' ,/  , '  .SET  DSKPRI  O', //.RUN  PHASE'/ 

2  '*' ,I/,'*',F/,'*' ,F/,'*' ,F/,'*Y' ,/ 

3  '*' ,1/,'*' ,F/,'*',F/,'*' ,F/,'*N') 

IF(NLINE.EQ.l)  GOTO  1100 

WRITE( 2 1 , 1035) 

1035  FORMATC' .RUN  BATCEN'/, 

1  '.SUBMIT  RDN . CTL/ RESTARTABLE/ UNIQUE : 0 ' ) 

1100  CALL  CLOSE( 21) 

STOP 

END 
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Q* ************** *************** *************** ******* *********** ******** 

c 

C  PHASE-CANCELLATION  STUDIES:  BIOENGINEERING 
C 

C  CALCULATE  HISTOCRAM  OF  AVERAGE  RECEIVED  PRESSURE  DATA  FOR  A 
C  RANDOM  PLANAR  DISTRIBUTION  OF  SCATTERERS 
C 

C* ******************************** ************************************** 


DIMENS ION  AMP( 4 1 , 4 1 ) , THETA( 4 1 , 4 1 ) , PIC( 4 1 , 4 1 ) , AVGDAT( 1 / 300 ) 

COMPLEX  PRS( 41 , 4 1 ) ,  C 

C  INPUT  VARIABLES  AND  INTIALIZE  ULTRASONIC  PARAMETERS 
TYPE  100 

100  FORMAT( 'OPLEASE  ENTER  RANGE, APERT, NUM:  ',$) 

ACCEPT  105, RANGE, APERT, NUM 
105  FORMAT(F,F,I) 

FREQ  ■  5.00 

M  =  41 

DIM  »  23.0 

SCAT  =*  1.0 

VELCTY  -  1500.00 

WAVEL  =  ( VELCTY/ ( FREQ* 1.0E03)) 

XK  =  (2.0  *  3. 14 159) /(WAVEL/ 10.0) 

CENTER  =■  INT(M/2 .0)  +  1.0 

TRANR  -  APERT/ 2.0 

ASF  =  (APERT/ (DIM-1 .00) ) 

DO  1005  I  -1, M 
DO  1000  J  -1,M 
1000  PIC(I.J)  » 

1005  CONTINUE 

C  DETERMINE  PARTICLE  LOCATION  IN  UNITS  OF  BLOCKS  AND  SEE  IF  WITHIN  BEAM 
C  THEN  SEARCH  TRANSDUCER  COORDINATES  -  CHECK  IT  WITHIN  BEAM  ALSO  - 
C  CALCULATE  THE  DISTANCE  BETWEEN  THE  TWO  SETS  OF  COORDINATES  AND  CONVERT 
C  CENTIMETERS  -  USE  THIS  VALUE  IN  THE  GEOMETRY 

DO  1045  I  =  1 ,NUM 
1010  SEED1  =  SEED1  +  .131 
SEED2  =  SEED2  +  .1331 
XCOOR  =  RAN(SEEDl)  *  FLOAT(M) 

YCOOR  -  RAN( SEED2)  *  FLOAT(M) 

TEMPI  =  ASF  *  SQRT((XC00R-CENTER)**2  +  (YCOOR-CENTER)**2) 

IF ( TEMP l.GT. TRANR)  COTO  1010 
DO  1020  II  =*  1  ,M 
DO  1015  J1  =  1 ,M 

T=ASF*( ( I1-CENTER)**2  +  (J 1-CENTER) **2) 

IF(T.GT. TRANR)  GOTO  1015 
PIC(II,J1)  -  'X' 

DIST»ASF*SQRT( ( XC00R-I1 )**2  +  (YC00R-J1)**2) 

HYP0“SQRT( DIST**2  +  RANGE** 2) 

PRS( 1 1 , J 1 ) - ( S  CAT/ HYPO) *CEXP( ( 0 . 0 , 1 . 0 ) *XK* ( RANCE+HYPO ) ) 

1015  CONTINUE 
1020  CONTINUE 
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C  LOCATE  CENTER  VALUE  OF  THE  PRESSURE  ARRAY  AND  ASSICN  TO  C  FOR 
C  NORMALIZATION  OF  THE  PRESSURE  VALUES  -  CALCULATE  THE  AMPLITUDE  AND 
C  PHASE  VALUES  AT  EACH  LOCATION  ON  THE  TRANSDUCER  SURFACE  AREA 

C  =  PRS( IF IX( CENTER) , IF IX( CENTER) ) 

DO  1030  II  =  10,32 

DO  1025  JJ  =  10,32 

IF(PIC(II,JJ).EQ./-')  GOTO  1025 
PRS(II.JJ)  =  PRS(II,JJ)/C 

THETA( 1 1 , J J )  =  5  7 . 29  580*  ATAN2  (  AIMAG(  PRS(  II,  JJ)  )  , 

1  REAL( PRS( II, JJ) ) ) 

AMP(II.JJ)  =  COSD( THETA( I I , J J ) ) 

1025  CONTINUE 
1030  CONTINUE 

K  =  0 

TEMP  =0.0 
DO  1040  II  =  10,32 

DO  1035  JJ  =  10,32 

IF(PIC(II,JJ).EQ.'-')  GOTO  1035 
TEMP  =  TEMP  +  AMP(II.JJ) 


K  =  K  +  1 

1035 

CONTINUE 

1040 

CONTINUE 

AVGDAT(I)  =  TEMP/ FLOAT( K) 

1045 

CONTINUE 

C  END  NUMBER  LOOP  AND  CALCULATE  MEAN  AND  VARIANCE  OF  AVGDAT  ARRAY 
K  =  0 

TEMP  =0.0 
DO  1050  I  =1 ,NUM 

TEMP  =  TEMP  +  AVGDAT(I) 

K  =  K  +  1 

1050  CONTINUE 

AMEAN  =  TEMP/ FLO AT(K) 

TEMP  =0.0 
DO  1055  I  =1  ,NUM 

TEMP  =  TEMP  +  (AVGDAT(I)-AMEAN)**2 

1055  CONTINUE 

AVAR  =  TEMP/ (FLOAT(K)-i .0) 

C  WRITE  RESULTS  TO  DISK  IN  FILE  SUM. DAT 
CALL  0PEN(21 , ' FILE' , ' SUM. DAT' ) 

WRITE( 2 1,1060)  AMEAN, AVAR 

1060  FORMAT( 5X , ' MEAN*  '.F11.7,'  VARIANCE  =  ' ,F11 .7,//) 

WRITE( 2 1,1061)  RANCE, APERT, NUM 

1061  FORMAT(5X,' RANGE  =  ',F7.2,'  APERTURE  =  ',F6.2,'  NUM  *  ',13,//) 
DO  1070  I  =1 ,NUM 

WRITE(21 , 1065)  I,AVCDAT(I) 

1065  FORMAT( 5X,I3,3X,F11.5) 

1070  CONTINUE 

CALL  CL0SE(21) 

STOP 

END 
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